{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "import scipy.stats.stats as scss\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# 设置随机数种子，使运行结果可复现\n",
    "np.random.seed(2046)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.894992</td>\n",
       "      <td>0.297397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.596078</td>\n",
       "      <td>-0.888134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.435564</td>\n",
       "      <td>-1.979664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1.153073</td>\n",
       "      <td>1.115046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.340793</td>\n",
       "      <td>0.064129</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1.621597</td>\n",
       "      <td>-0.930433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2.139849</td>\n",
       "      <td>2.060672</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.488654</td>\n",
       "      <td>1.067725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2.049478</td>\n",
       "      <td>-0.115373</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.564523</td>\n",
       "      <td>0.206383</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.270165</td>\n",
       "      <td>-1.018755</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.596075</td>\n",
       "      <td>-2.018875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1.291240</td>\n",
       "      <td>1.134136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.251871</td>\n",
       "      <td>0.026758</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1.110560</td>\n",
       "      <td>-1.144212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2.805694</td>\n",
       "      <td>2.271965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.465100</td>\n",
       "      <td>1.058364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2.258865</td>\n",
       "      <td>-0.081940</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    x1  x2        x3         y\n",
       "0    0   0  0.894992  0.297397\n",
       "1    0   1  0.596078 -0.888134\n",
       "2    0   2  0.435564 -1.979664\n",
       "3    1   0  1.153073  1.115046\n",
       "4    1   1  1.340793  0.064129\n",
       "5    1   2  1.621597 -0.930433\n",
       "6    2   0  2.139849  2.060672\n",
       "7    2   1  2.488654  1.067725\n",
       "8    2   2  2.049478 -0.115373\n",
       "9    0   0  0.564523  0.206383\n",
       "10   0   1  0.270165 -1.018755\n",
       "11   0   2  0.596075 -2.018875\n",
       "12   1   0  1.291240  1.134136\n",
       "13   1   1  1.251871  0.026758\n",
       "14   1   2  1.110560 -1.144212\n",
       "15   2   0  2.805694  2.271965\n",
       "16   2   1  2.465100  1.058364\n",
       "17   2   2  2.258865 -0.081940"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 生成模型数据，其中x1，x2为不相关的特征；x1，x3强相关\n",
    "n = 2\n",
    "data = []\n",
    "for i in range(0, 3):\n",
    "    for j in range(0, 3):\n",
    "        data.append({'x1': i, 'x2': j})\n",
    "data = pd.DataFrame(data * n)\n",
    "# 生成强相关自特征\n",
    "data['x3'] = data['x1'] + np.random.random(9 * n)\n",
    "# 生成被预测值\n",
    "error = 0.1 * np.random.random(9 * n)\n",
    "data['y'] = 0.7 * data['x1'] - 1.1 * data['x2'] + 0.3 * data['x3'] + error\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABaUAAAH3CAYAAACrTTPQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAA9hAAAPYQGoP6dpAABTeElEQVR4nO3de3hdZZ03/O9KSxvANlqwNohCwcMYqkjVjh3wOINWxjozoo7PWB1HQcXD4PHxONZ6GHV0RnR0UFFEAcfjqG8VKo+CiIhWKSqloqIBBYodCIYihLbJ/f6RJjYkadNmZ+0cPp/rypVmrXvv9Uv26v7t9V2nqpQSAAAAAACoQ0uzCwAAAAAAYOYQSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALWZ3ewC9qSqqirJIUm2NrsWAKateUluLKWUZhcyVenXANRAvx4n/RqAmuyxZ0/6UDr9DfP6ZhcBwLR3aJIbml3EFKZfA1AH/Xp89GsA6rLbnj0VQumtSfK73/0u8+fPb3YtAEwzt912W+53v/sljhgaL/0agAmjXzeMfg3AhBprz54KoXSSZP78+ZomAExy+jUATH76NQDN5kaHAAAAAADURigNAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURigNAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADUZnazC6jDth19Ofuya3Nd1x05bMEBee7ywzNntjye6ce6zkxgPQeAqaG3r2R9Z1e2bO3JwnmtWbZ4QWa1VM0uiwbw2gIwXtM+lH73eZtyxiWd6St/mvau836ekx+zOG88oaN5hUGDWdeZCaznADA1rNu4OWvWbsrm7p7Bae1trVm9siMrlrQ3sTLGy2sLQCNM60PL3n3epnzsu0PDiyTpK8nHvtuZd5+3qTmFQYNZ15kJrOcAMDWs27g5p5yzYUhomSQ3dffklHM2ZN3GzU2qjPHy2gLQKNM2lN62oy9nXNK52zFnXNKZbTv6aqoIJoZ1nZnAeg4AU0NvX8matZtSRpg3MG3N2k3pvfteZiY9ry0AjTRtQ+mzL7t22NF0d9dX+sfBVGZdZyawngPA1LC+s2vYUbS7Kkk2d/dkfWdXfUXREF5bABpp2obS13Xd0dBxMFlZ15kJrOcAMDVs2Tp6aLkv45g8vLYANNK0DaUPW3BAQ8fBZGVdZyawngPA1LBwXmtDxzF5eG0BaKRpG0o/d/nhaal2P6al6h8HU5l1nZnAeg4AU8OyxQvS3taa0dp2laS9rTXLFi+osywawGsLQCNN21B6zuyWnPyYxbsdc/JjFmfO7Gn7J2CGsK4zE1jPAWBqmNVSZfXKjiQZFl4O/Lx6ZUdm7WlvM5OO1xaARprWW+9vPKEjL37s4mFH17VUyYsfuzhvPKGjOYVBg1nXmQms5wAwNaxY0p7TVy3Norahl3FY1Naa01ctzYol7U2qjPHy2gLQKFUppdk17FZVVfOTdHd3d2f+/Pn79BzbdvTl7MuuzXVdd+SwBQfkucsPdzQd05J1nZmg0ev5bbfdlra2tiRpK6Xc1rBCZ5hG9GsAppfevpL1nV3ZsrUnC+f1X9ZhX4+i1a8bo1H9upGvLQDTy1h79owIpQFgNDZyG0O/BmAi6deNoV8DMNHG2rMdQgkAAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURihNw/T19eUJT3hCHv/4x+eSSy6pbbnr16/PTTfdtNsxn/3sZ3PDDTfUVNGf9PX15aqrrsonP/nJnH766Xv12D/84Q/58Ic/nI985CO7HXfllVemq6trPGUCQBK9fCR6OQDNoi8Ppy/DNFJKmdRfSeYnKd3d3YXJb7/99itJynnnnVfL8s4888yy3377lSc84Qmlt7d3xDHbt28vBx98cGltbS2vfvWrG15DX19fueWWW8oVV1xRvvzlL5f3v//95aUvfWl54hOfWObPn1+SDH595StfGfV5zj///PIP//APZceOHaWUUq655pqSpMydO3dwzLnnnltOPPHE8sMf/nBw2nOe85yy3377lZNPPrnhvxvMBN3d3QP/R+eXSdD3puqXfj196OV6OUxG+rV+PVPpy/oyTDVj7dmz6wi+mTnmzp2b7du3Z86cOcPmPeABD8ivf/3rMT/Xi170onzsYx/b7ZjDDjssO3bsyEUXXZR3vetd+Zd/+ZdhY771rW/l5ptvTpIcccQRY15+d3d3Xvayl6W3tze9vb3Ztm1benp60tPTk9tvvz1bt27NH/7wh3R1dWXHjh0jPkdVVTn44IOzcOHCzJs3L5/61KfymMc8JgcddNCQcTt27Mhb3vKWXH755bnnPe+Zj3zkI9lvv/2S9P9NB3z0ox/NJZdckhe+8IWD0zZu3Jjt27fngQ984Jh/NwAYjV4+lF4OQDPpy0PpyzB9CKVpqFmzZiVJZs8evmq1trYmSdrb23PAAQeM+hw333xzuru7R2y6d/fEJz4xr3zlK/OBD3wg73jHO/K0pz0tRx999JAxZ511VpJk6dKlOeWUU8b6q6StrS2dnZ35/ve/P+LvcsABB6Snpyc7duzI8ccfn0c84hFpb2/Pfe973xx66KE59NBDc5/73GfEv8XdzZ49O9/85jezfPny/Nd//Vee8IQn5LjjjkuSwcb5wx/+MJdcckmOPfbYPOUpT0mS9PT05KqrrkqSPPOZzxzz7wYAo9HL9XIAJg99WV+GaWt3h1FPhq84vWhKude97lWSlEsuuWTYvEc84hElSVm7du1un+PUU08tScqrXvWqMS3ztttuK/e5z31KkvIP//APQ+bdfPPNZe7cuSVJueCCC8b+i+z0q1/9qlx22WXl5z//ebnxxhtLd3f34Kk/119/fVm4cGFJUg4++OBy+eWX7/Xz393FF19cjj766HLppZeWW265pSQp9773vUsppbz5zW8uVVWVr33ta4PjL7roopKkHHXUUeNe9mR1zTXXlFWrVpWFCxeWuXPnlqOPPrp85jOfmbDlffGLXyx/8Rd/UQ444IAyf/788nd/93fl6quv3u1jtm3bVt797neXhzzkIWXOnDmlvb29/PM//3P5wx/+MGF10jhOB9avGUovH5/p3Mtvuumm8pKXvKQceuihZc6cOeXBD35w+cAHPjDq6d3j8fvf/76sWbOmLFu2rMybN68ceOCB5dGPfnT58Ic/XHp6esb8PO94xztKknLQQQeVzZs3N7xO6qNf69czlb48PtO5L++rybyN/bjHPW7I5Vl293XYYYdNWM2Mj8t30BRVVY06b2AP71gN7L3ck3nz5uU973lPbr311vzzP//zkHmf+MQnctddd+Vxj3tcjj/++L1aftJ/OtQDHvCAYdO7urry9Kc/PVu2bMmhhx6aM888M0uXLh0y5sILL8xll12W173udXvcI/3Nb34zGzduzNy5c/OSl7wkmzZtyoYNG5Ikd911Vz7xiU/k8MMPzxve8IZ0dnbmve99b1760pfm4osvTpI8+clP3uvfbSq44oor8vjHPz633XZbkv7166c//Wme97zn5ZprrsmaNWsaurx/+Zd/yTvf+c7Bn6uqyle+8pVceOGFueiii3LMMccMe8xdd92Vpz71qfnWt741+JjNmzfnQx/6UC688MJceumlmT9/fkPrBJhIevmf6OV/ct111+W4447L9ddfn6R/PfnFL36RV73qVbniiivy6U9/umHLuvjii/P0pz998CZTs2bNSm9vb37wgx/kBz/4Qc4888x8/etfT3t7+26f54orrsjb3/72JMlHPvKRLFq0qGE1AtRFX/4TfXn8Jvs29pw5c4ZcXmUkvb292bFjR1paWhpaK02wu8R6MnzFntxJ61vf+laZO3duude97lXuc5/7lEMOOaS0tLQM7nk85JBDysKFC8v8+fPLxRdfXJYvX75Xe3Hf8pa3DJv3xS9+sSxYsKDc7373K4sXLy5HHnnkbr8OOOCAwaNjBqYdccQR5f73v3+5z33uU5YtW7bXv/evfvWr8uAHP7gkKSeeeGLp6uoaMr+np6e89rWvLVVVlSTl/ve/f/nWt7612+d82cteNua9gQNfv/vd7wb/posWLSoPfvCDh30dfvjh5ZBDDtnr33Ey6OrqGtw7v3jx4nLxxReX3t7esn79+nLYYYeVlpaW8v3vf79hyzv33HMH/7bPec5zyo033lh6enrKBz7wgVJVVXnoQx9atm3bNuxxJ598cklSWlpayr/927+V22+/vdx8881l1apVJUl5yUte0rAamRiOvNKvZzK9XC8fi23btpWOjo7B1+GrX/1q2bFjR7n66qvLwx/+8JKkfO5zn2vIsn77298OHhX4kIc8pFx44YWlt7e33HrrreW0004bPDrvmGOOGTyybiQ9PT3lqKOOKknKM5/5zIbURnPp1/r1TKAv68sTaapsY+/JCSecUJKU1atXN6xWGsuR0ky4/fffPw94wANywAEHZNasWbnlllvS19eXJFm0aFHucY97ZMeOHenp6Ulra+te78Uaaa/vjh070tXVNXjkzFjdcsstueWWW4ZNP/jgg8f8HDfeeGPOOOOMvO9970tfX19OP/30vOQlLxky5tJLL81JJ52Uq6++Okly7LHH5i1veUv+8i//crfP/eQnPzkHHXRQ5s6dm7lz5+ass87Kxo0bk/TvpX7Xu96VUkp6e3uzffv23HHHHbnjjjuyfv36JMlNN92Um266acTn3tNexsnqX//1X/P73/8++++/f775zW8O3mTiUY96VM4+++w89rGPzamnnjr4NxiPnp6evOENb0iSPOlJT8rZZ589eETCK1/5ylx55ZU588wzc8YZZ+SlL33p4OOuvPLKfPKTn0ySvP3tb8/rXve6JMmBBx6YM888Mz/60Y/ysY99LC972cuyZMmScdcJ0Gh6uV4+Fh//+MezadOmwaObHvOYxyRJHvzgB+fLX/5yHvzgB+eVr3xlnvGMZ+z1UXt39973vje33npr7nnPe+bb3/724NHQ97znPXPqqadm3rx5eeELX5grrrgiX/3qV3PiiSeO+DxvfvObc9VVV+U+97lPTj/99HHVBFAXfVlfnkhTYRt7Ty644IKcd955ud/97pfXv/71466TJttdYj0ZvmJP7pSxevXqwb1gI13v6thjjy1JyqGHHjriHseBr4GjY0ba63XnnXeWW265pfT09Ix6dMzAXuDR9pr19fWVu+66q2zdurVs2bJlj79XX1/f4JE2A1+PfOQjy3Oe85whX0972tMG99z+2Z/9WTn//PP3+Nwj+drXvlaqqir3vOc9S5Jyn/vcp2zdurW8853vLH19fYPjPvShD5Uk5QlPeMKw5zjiiCNKknLdddcN28s8FWzfvn3w93/FK14x4phHPepRg7/jeH3pS18afG1Hum7ZVVddVZKUJz7xiUOmv/zlLx+83tldd9017HEf+chHSpLy9re/fdw1MnEceaVf8yd6uV4+kiVLlpQkZeXKlSPOf+Yzn1mSlO9+97vjXtbAEVyvec1rRpy/ffv2wfXrta997YhjLr744sEjC7/61a+OuyYmB/1av56J9GV9uVGmyjb27uzYsWPwM8l///d/j7tGJo4jpandD3/4wzGNG7gW4b5obW0dvMPwvqqqKnPmzMmcOXNyj3vcY0zjOzo6Bu/AmyQ//vGP8+Mf/3jY2AMOOCBvfetb8+pXv3rM1+va1Ve+8pX8n//zf9LS0pKPf/zjedaznpW+vr685S1vyQc/+MHceuutef/7358k+dznPpckw+6EnGRwL/dBBx2UAw88cK/raLYrr7wyf/jDH5Ikz3rWs0Ycc8IJJ+RHP/pR1q1blxe96EXjWt53v/vdJMmRRx457LplSdLR0ZHDDjssl1xySf74xz8O/k0HHvc3f/M3I17T7IQTTkiSnH/++fmXf/mXcdUIUAe9XC+/u1tvvXXwyLLd9eQvfvGLOf/88wePot4XW7duze9///skyeMf//gRx8yePTuLFi3Krbfemu7u7hGf4/nPf376+vryvOc9L3/zN3+zz/UANJu+rC83ylTZxt6dM844Ixs3bsxxxx2XZz/72eOqj8nBVcFpiG3btuWSSy4ZMu32228fcezatWt3u/f+1FNPraPkvfKGN7whX/ziF/PjH/84N9xwQ26//fbccsstOfnkkwfHPOUpT8lVV12V17/+9dlvv/1y/fXX55nPfObgm/Hu9Pb25m1ve1ue+cxnZvv27TnrrLMGT0fasWNH3vCGN+TQQw/Nv//7v+e0005LZ2dnLrvssiTJL37xi2HP1d3dnQMPPHBKNsskueaaa5L0n172qEc9asQxD3vYw5Jk8DSuRixv+fLlo4552MMelu3bt+fXv/71mB93+OGHZ/78+Q2pEWCi6eV6+UgGel0yer9rVE/eNWzYf//9Rx138803J0kOOeSQYfNe85rXpLOzM3PmzElHR0de+9rX5tRTT81HP/rRvT41HaCZ9GV9uZGmyjb2aLZt25Z3vOMdSTL4nalPKE1DfO9738sf//jHwZ9/8pOf5KEPfWhuuOGGJlbVOEuXLs0znvGMPOIRj8ghhxyS73znOznmmGNyxhln5JBDDsnnP//5nHfeeTn88MOT9DetF7zgBfnSl76Uv/qrv8qZZ5456nNv3Lgxj3rUo7JmzZrMnj07n/vc57Jq1arceeedSfrffBctWpQvfelLmTNnTj7zmc/kfe9738Dpd7nyyiuHPN8NN9yQUkruf//7T8wfowYD1yY77LDDRr1e18KFC5Mk1157bcOW96AHPWjUMXdfXk9PT+644449Pu7e9753br311sG7GwNMVnq5Xj6SgR653377Df5t7q5RPbm1tTWLFi1KklxxxRUjjrn66qvzv//7v0mSFStWDJl36aWX5hOf+ESSZPv27fn85z+fq666Kt/85jdzyimn5Mgjj8w555wzrhoB6qIv68uNNBW2sXfnrLPOyo033phHP/rRo55NxdQjlKYhvvCFLwz5+e1vf3uuvfbaPOMZz8i2bduaVFXjXX755Xnyk5+cpz71qbnhhhvy8pe/PD//+c+Hnf4ya9asfOMb38gLXvCCbN++PS984Qvz9re/fcTnPOqoo3L88cfnsMMOy0UXXZRnPvOZ2bZtW37961/nxBNPzF//9V8nSf78z/88n/3sZ3PuuefmrLPOSpLc5z73yfXXX58tW7YMPt/vfve7JMkRRxwxpt9p9erVOemkk3LSSSflfe97397+SXbrSU960uDpYGP5GtgrftdddyVJ7nWve4363APzNm/ePO4692V5A4+ps06AiaSX6+UjGeh38+fPH/Umho3sdc94xjOSJB/+8Idz6623DpnX29ubV7ziFUmShz/84Xn0ox89OG/gSMBSShYuXJgf/vCH2bBhQ84///xcffXVueiiizJ37tw897nPzWc+85lx1wkw0fRlfXkk03kbezS9vb35t3/7tyQZvHki04NQmnHr6enJl770pSQZvPvvv/3bv+Xggw/OD37wg7zpTW9qZnmj6u3tTU9PT2655ZZRm/qWLVvy9a9/Pe985ztz9NFH55GPfGQuuOCCHHvssfnxj3+c//zP/8z8+fMHT+e58cYb86tf/So/+clPsn79+vz93//94Okqq1evHtxY2lVVVXn3u9+dn/3sZ4NjL7nkkjzucY/Lb37zm/zHf/zH4NgTTzwx5557bu6888487nGPy9/+7d8myeBpRknyq1/9KkkG76S7J1/+8pfzyU9+Mp/85CfzjW98Y2x/vDHatm1b7rrrrjF/bd++Pcmf7gp9wAEHjPrcA9dwHtjbPR77srxdN8zrqhNgoujlevlo6u7Jb3/72/OQhzwk1113XZ70pCfl+9//frq7u/P9738/T3va0/Ktb30rSYbdq+Giiy7K5ZdfniT54Ac/OOzU5Mc//vH5+Mc/niT5v//3/6a3t3fctQJMFH1ZXx7NdN7GHs0Xv/jF/PrXv05HR0ee9rSnjbs2Jg83OmTcPvvZz+aWW27J0qVLc9111+WWW27J4Ycfns9//vP5q7/6q3zgAx9oyE1m/umf/ilf+9rXMnfu3FGP1EkyeNObNWvWZM2aNWN67ksuuSTHHXfcsOnf/OY387znPW/ItNmzZ+e2227LiSeemK1bt+b2228f85v2hz70ofzhD3/ImWeeOeR3uPnmmzN79uz88Y9/zJw5c3LBBRck6W96973vfYc8x8tf/vLccccdefKTn5wtW7bkYx/7WC688MLBv/GPfvSjJMlDH/rQMdU0kb7zne/s0+MGriM50s0DBwx8OBu4hMZ47Mvydr3WZV11AkwUvVwvH03dPfle97pXLrvssrzxjW/Mpz/96Rx77LHDxhx77LF5+tOfPmTahRdemCQ5+OCDc+KJJ4743E95ylPS2tqa3//+9/nhD3+Yv/iLvxh3vQATQV/Wl0cznbexR/Of//mfSZKXvOQlqapq3LUxeQilGZe+vr7B0yj+8R//ccjpM0984hNz0kkn5Ywzzsj3vve9wekrV67cp2XNmjUrs2bNyn777Tf475EMXGcw6b/T62h6enpy1113paenZ/DN8O5WrlyZOXPmZMeOHTnooIOyYMGCzJ8/PwceeGDud7/7Ze7cuYMN7tZbb82KFSvygAc8IPPmzcu8efNyj3vcY/BmCOeff34+/elP5zOf+UwOPvjg/Pu///vgch72sIcN3m1+V1/4wheGnbaVJD/72c/y0Ic+NDfeeGOS5Otf/3o++MEPJknWr1+fJDnmmGNG/d0nu4HTeHZ3vbSBU3rvvle8ruXNmjUr8+bNy9atW3PDDTfk0EMPnfA6ASaCXq6X785AjxyocySN7nVtbW35r//6r5x22mn5xS9+kTvuuCPnnntu/vM//zMtLS057bTThj1mYJ055phjhtwwcVf77bdfFixYkBtvvDG//e1vhdLApKQv68sTYSpsY4/kl7/8Zb7//e9n9uzZefaznz3uuphchNKMy7nnnptf/OIXucc97jGsYSbJe9/73jzykY/Mi170osHTVtrb23d7CsfNN988uCd2V5/4xCcGb14zmltuuWVIOLhmzZo85znP2ZtfaYh73vOe+c1vfpOFCxeOuoGT9F/X8NZbb83rX//6US+6/6xnPSs9PT253/3ul3e+851D5j3wgQ/MIYcckv333z8333xzfvnLX+be9753Ojo6hvxuGzduzIMe9KDBPbSHHHJIli1blvXr1+cHP/hBDjvssGzYsCELFizI0UcfPabfcePGjWMaV6eBG0hce+21KaWMuDd04APGPe5xj4Ytr7Ozc9QxIy3v/ve/f6666qp0dnbmz//8z8f8OIDJRC/vp5ePbKBH3nnnnbnpppsGb0S4q4nqdXPmzMlDH/rQdHd359xzz02SvOAFL8gjH/nIYWPnz5+fpP/13p2tW7cmSQ488MCG1grQKPpyP325sabKNvbdfepTn0rSf3Pje9/73uOui8nFNaXZZ1u3bs2b3/zmJMnJJ5+ctra2YWPuda975UUvetGQaR//+MdzzTXXjPr1/Oc/f59r+tCHPpSenp7BN7VXvOIVueaaa/b5+ZLkvve9b/bbb7/85Cc/GbKHeCy+9KUvDTb/qqry3//933nf+9437G63l1xySTZs2JBLL700D3nIQ5L0n6Lyne98Z/BroBHf/YYPJ510UpLkAx/4QL7whS+kr68vT3rSk0bdMz0VHHXUUZk7d256enpy1VVXjTjmBz/4QZI05A7Ij3jEI5IkP/7xj0cdM9Ly9vS4X//617n55psbVidAo+nlezbTe/nBBx882MNG63eN7Mkjee9735uurq4sWLAg73nPe0Ycc9RRRyVJrr766lGfZ9OmTYOh9MBrBDCZ6Mt7NtP78r6aKtvYu+rt7R28OfFzn/vccdfEJFRKmdRfSeYnKd3d3YXJ5V3veldJUubNm1e2bNlSSinloIMOKknKRRddNGz8scceW5KUtWvX7vZ5Tz311JKkrF69eq/q+cUvflFaW1tLkvKpT32qnHjiiSVJedCDHlRuvPHGvXquu7v22mvLwoULy6JFi8oXv/jFYfOPPvroYb/3WWedVVpaWsrixYvL+vXrx7Scyy+/vCQpD3/4w0tfX9/g9Jtvvrnc4x73KPvvv3/ZvHnzkMfccccd5aCDDipVVZVDDz20JClf//rX9+0XnUT+8i//siQpb37zm0ec/+d//uclSXnb29427mX99re/LVVVlaqqyi9/+cth83/+85+XJCVJufbaawenn3322SVJOeyww8qOHTuGPe70008vScrhhx8+7hqZON3d3QOv7/wyCfreVP3Sr6cmvfxP9PLRvfCFLyxJynOe85wR5//93/99SVKe//znN3zZ119/fTnggANKkvLxj3981HE33nhjmT179m7/dgN1/tmf/VnD62Ti6df69UygL/+Jvtx4U2Ebe1cXX3xxSVJmz55dtm7dOu6aqM9Ye/bU3c1D073qVa/KMccckzVr1jT9NIpbb701f/M3f5Oenp78xV/8RZ73vOflE5/4RA477LD88pe/zHHHHZdNmzbt03Nv2bIlT3nKU7Jly5bceeed6erqGtPj9t9//7S1taWzszPHHnvskDv8juZ1r3tdkuSaa67JX//1X+e9731vLrvssrz61a/O7bffnlNOOWXYabP7779/XvziF6eUkuuvvz5HHHFEVqxYMebfb/Xq1TnppJNy0kkn5X3ve9+YHzfRBvbmf/jDH871118/ZN6FF16YH/7wh0mSpz71qeNe1v3ud7888YlPTCklb3jDG4bN/9d//dckyZIlS3LYYYcNTv+7v/u7zJ8/P9ddd13+67/+a8hjenp6Bl/zRtQIMBH08t3Ty/sN9OTPf/7z2bBhw5B5V199db7yla8kmZh+9+pXvzp33HFHli9fPnjk2kja29vz0pe+NEmyatWqnHXWWdm6dWtKKdm0aVOe9axn5fOf/3yS/r8XwGSkL++evjw+U2Ebe1fnnXdekuRRj3qUy2FOV7tLrCfDV+zJndSuu+660tvbO/hzI/fivulNbxpTDTfccEN56EMfWpKUQw45pFx33XWD8y6//PJyr3vdqyQpBxxwQPnIRz4ypN49uf7668tRRx1VkpQjjjii/OIXvxhx3Eh7cUvp3/u7bNmywT2AT3/608ttt9026vKuuOKK8t73vrc86UlPGjwqaNev0047rdx1113DHnf++ecPjjn55JPH/PuVUgZ/vyTlcY973F49diL19PSUxYsXlyTlwQ9+cPne975X7rzzzvI///M/ZcGCBSVJeexjHzvsce94xzvK3Llzy9y5c/dqeRdccMHg3+EFL3hBufHGG8uWLVvKq171qsHpZ5555rDHvelNbxrce3vaaaeV2267rVx99dWDe6Fnz5496nrD5ODIK/16ptPL++nluzfw2i9atKh84xvfKD09PeXCCy8shx9+eElSjjzyyLJt27Yhj/nMZz4z2JNHOwpqd7797W8P9tKf/vSnexz/xz/+sTzjGc8Y8vduaWkZ8vOpp56613UwOejX+vVMoS/305cbb6psYw8YWAfHut4yeYy1Zze9Ke7pS9OcWnbXMB/96EfvVcN83etet8flffnLXy4HH3xwSVIOPfTQctVVVw0b87Of/awsWrRo8E3v6KOPLl/4whf22Dh/9KMflfve974lSVm6dOng6VMjGWg6I/3ed911V3nuc587uPxjjjmm/P73v9/j73bXXXeV17/+9cOa5sKFC8tZZ501OG79+vWDDSRJqapqt2/so9U+2RpmKUM/8Nz9a8GCBSN+gFm9evXgmL311re+dcRlJSlPe9rTRlxnenp6ypOe9KRRH/f+979/n3536mMjV79mKL38omHz9PL+EOD+97//iL2utbW1fO973xv2mE996lODYzo7O/dqedu2bSsPechDSpLy2te+dq8eu3bt2nLiiSeWI444orS2tpYDDzywLF++vHz605/eq+dhctGv9euZSl++aNg8fXnfTYVt7FJK+d3vfjc47v/9v/+318uluYTSNMU973nPkqR8+9vfHjbvUY961KhvRiN9vfKVrxx1OZdeeml54hOfODj2uOOOK9dff/2o46+55prBvWwDX/e9733La17zmvKd73yn9PT0DBn/0Y9+tOy3334lSXn84x8/bP3bsWNHueOOO8ott9xSvve975W5c+eWJOXSSy8dtYY3vOENJUlpa2srl19++ajjSuk/0uetb31rmTVrVtlvv/3KueeeWy6++OKyatWqMmvWrHL++eeXvr6+ctppp5U5c+aUJOXFL35xecc73jHYNN/97nfvdhkDJnPDLKWU3/zmN+Vv//ZvhxzpdPzxx5df/epXE7K8L33pS+WBD3zg4LIOPPDA8qY3vWnY0V+72r59e3nnO9855IPLoYceWs4555wJqZHGspGrXzOUXq6Xj2bLli3l+c9//mC9ScqyZcvGfF1PGA/9Wr+eqfRlfbnRpsI2NlObUJqmOPDAA0uScv755w+b9/CHP7wkKe3t7eXII48c9autra0kKS972cuGPP7aa68tp5122uDzJP2nC73nPe8Z8SZzd7dt27bB007u3pz/4z/+Y8jYW2+9tZxwwgnlaU97WrnzzjuHPdf3vve9Yc9RVdUe986uWbOm/PCHPxx1/u9+97uyZs2acsghh5Qk5YEPfOCwI4+uu+66smHDhrJ8+fLB5a5evXrwpg3/9E//VEtjqduWLVvK+vXrx31DjbHo6+srV199dbn88stHfP1H09PTU6644opy1VVXDbmJBpObjVz9mqH0cr18T/7whz+UH/3oR/t0SQ7YV/q1fj1T6cv68kSZCtvYTE1CaZpi4K7nX/va14bNGzgFc6ynFr3oRS8aMv0lL3nJYCNoaWkpz3ve8/ZpY+jqq68u//iP/zh4F+G//Mu/HHHc9u3by/bt20d9nkc+8pGD9dz//vcvH/zgB/e6lgGnnXZaecQjHlGqqipJ/2kz73rXu8odd9wxZNxNN91U/vZv/3Zw3OLFi8sFF1wwZExvb2/553/+58HaDj744PK73/1un2uD6c5Grn7NUHr5vtHLYWLp1/r1TKUv7xt9GZpHKE3ttm/fPvgm/bnPfW7Y/COOOGJMDXN3/uqv/qqcfPLJ5eqrrx5PqaWUUv73f/+3/Ou//mu58sor9+nx3/72t8tXv/rV8tvf/nbctfzqV78qBxxwQDnqqKPKBz7wgd3eqOEd73hHOfDAA8vq1avLH//4x1HHfe1rXyv3ve99yxVXXDHu+mA6s5GrX/Mnevm+08thYunX+vVMpC/vO30ZmmesPbsq/Y1p0qqqan6S7u7u7syfP7/Z5cCEuf7663PooYeOaeyWLVuycOHCPY7btm1b5syZM97SYFq77bbb0tbWliRtpZTbml3PVKVfg14OE0m/bgz9mplEX4bmGGvPbqmvJGB3xtosk4ypWSbRLAGgRno5AEwe+jJMbkJpAAAAAABqI5QGAAAAAKA2QmkAAAAAAGojlAYAAAAAoDZCaQAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2s5tdAAAAANNPb1/J+s6ubNnak4XzWrNs8YLMaqmaXRbAmHgPg4kllAYAAKCh1m3cnDVrN2Vzd8/gtPa21qxe2ZEVS9qbWBnAnnkPg4nn8h0AAAA0zLqNm3PKORuGhDlJclN3T045Z0PWbdzcpMoA9sx7GNRDKA0AAEBD9PaVrFm7KWWEeQPT1qzdlN6+kUYANJf3MKiPUBoAAICGWN/ZNezowl2VJJu7e7K+s6u+ogDGyHsY1EcoDQAAQENs2Tp6mLMv4wDq5D0M6iOUBgAAoCEOPnBuQ8cB1Ml7GNRHKA0AAEBjVA0eB1An72FQG6E0AAAADXHz7Xc1dBxAnbyHQX2E0gAAADTEwnmtDR0HUCfvYVAfoTQAAAANsWzxgrS3tY56ZnuVpL2tNcsWL6izLIAx8R4G9RFKAwAA0BCzWqqsXtmRZPglVwd+Xr2yI7NaXJAVmHy8h0F9hNIAAAA0zIol7Tl91dIsaht6evuittacvmppVixpb1JlAHvmPQzqMbvZBQAAADC9rFjSnuM7FmV9Z1e2bO3Jwnn9p7s7uhCYCryHwcQTSgMAANBws1qqLD/yoGaXAbBPvIfBxHL5DgAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2QmkAAAAAAGojlAYAAAAAoDZCaQAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2QmkAAAAAAGojlAYAAAAAoDZCaQAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2QmkAAAAAAGojlAYAAAAAoDZCaQAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2s5tdAAAAANNPb1/J+s6ubNnak4XzWrNs8YLMaqmaXRYAsItm9WuhNAAAAA21buPmrFm7KZu7ewantbe1ZvXKjqxY0t7EygCAAc3s1y7fAQAAQMOs27g5p5yzYcgGbpLc1N2TU87ZkHUbNzepMgBgQLP7tVAaAACAhujtK1mzdlPKCPMGpq1Zuym9fSONAADqMBn6tVAaAACAhljf2TXsiKtdlSSbu3uyvrOrvqIAgCEmQ78WSgMAANAQW7aOvoG7L+MAgMabDP1aKA0AAEBDLJzX2tBxAEDjTYZ+LZQGAACgIZYtXpD2ttZUo8yvkrS3tWbZ4gV1lgUA7GIy9GuhNAAAAA0xq6XK6pUdSTJsQ3fg59UrOzKrZbTNYABgok2Gfi2UBgAAoGFWLGnP6auWZlHb0FN+F7W15vRVS7NiSXuTKgMABjS7X8+e0GcHAABgxlmxpD3HdyzK+s6ubNnak4Xz+k8BdoQ0AEwezezXQmkAAAAablZLleVHHtTsMgCA3WhWv3b5DgAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2QmkAAAAAAGojlAYAAAAAoDZCaQAAAAAAajO72QUAAAAw/fT2lazv7MqWrT1ZOK81yxYvyKyWqtll0QBeWwDGSygNAABAQ63buDlr1m7K5u6ewWntba1ZvbIjK5a0N7EyxstrC0AjuHwHAAAADbNu4+accs6GIaFlktzU3ZNTztmQdRs3N6kyxstrC0CjCKUBAABoiN6+kjVrN6WMMG9g2pq1m9LbN9IIJjOvLQCNJJQGAACgIdZ3dg07inZXJcnm7p6s7+yqrygawmsLQCMJpQEAAGiILVtHDy33ZRyTh9cWgEYSSgMAANAQC+e1NnQck4fXFoBGEkoDAADQEMsWL0h7W2uqUeZXSdrbWrNs8YI6y6IBvLYANJJQGgAAgIaY1VJl9cqOJBkWXg78vHplR2a1jBZtMll5bQFoJKE0AAAADbNiSXtOX7U0i9qGXsZhUVtrTl+1NCuWtDepMsbLawtAo8xudgEAAABMLyuWtOf4jkVZ39mVLVt7snBe/2UdHEU79XltAWgEoTQAAAANN6ulyvIjD2p2GUwAry0A4+XyHQAAAAAA1EYoDQAAAABAbYTSAAAAAADUxjWlAQAAaLjevuJmeAAwyTWrXwulAQAAaKh1GzdnzdpN2dzdMzitva01q1d2ZMWS9iZWBgAMaGa/dvkOAAAAGmbdxs055ZwNQzZwk+Sm7p6ccs6GrNu4uUmVAQADmt2vhdIAAAA0RG9fyZq1m1JGmDcwbc3aTentG2kEAFCHydCvhdIAAAA0xPrOrmFHXO2qJNnc3ZP1nV31FQUADDEZ+rVQGgAAgIbYsnX0Ddx9GQcANN5k6NdCaQAAABpi4bzWho4DABpvMvRroTQAAAANsWzxgrS3taYaZX6VpL2tNcsWL6izLABgF5OhXwulAQAAaIhZLVVWr+xIkmEbugM/r17ZkVkto20GAwATbTL0a6E0AAAADbNiSXtOX7U0i9qGnvK7qK01p69amhVL2ptUGQAwoNn9evaEPjsAAAAzzool7Tm+Y1HWd3Zly9aeLJzXfwqwI6QBYPJoZr8WSgMAANBws1qqLD/yoGaXAQDsRrP6tct3AAAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFCb2c0uAAAAZpLevpL1nV3ZsrUnC+e1ZtniBZnVUjW7LAAAZqBmfTYVSgMAQE3WbdycNWs3ZXN3z+C09rbWrF7ZkRVL2ptYGQAAM00zP5u6fAcAANRg3cbNOeWcDUM+9CfJTd09OeWcDVm3cXOTKgMAYKZp9mdToTQAAEyw3r6SNWs3pYwwb2DamrWb0ts30ggAAGicyfDZ1OU7AGAGq6rqnklWJDksye1JrkxySSlFMgYNtL6za9hRKLsqSTZ392R9Z1eWH3lQfYUBU4aeDUCjTIbPpvsUSldVdVySZyQ5Mklvkp8n+UIp5YoxPPaBSb6ZpJRSjtyX5QMAu1dV1UOTrE7ymCStSa5O8oFSyud2GfO8JB9OcuDdHv6bqqpeVEq5qK56YbrbsnX0D/37Mg6YPvRsAOo2GT6b7lUoXVXVfkk+leT/3G3WyiT/t6qqLyd5aSnl5t08zZwkhycjHiEOAIxTVVVPSvK19PfcgdsmPzLJuVVVLSmlvKWqqr9O8skks0Z4iiOTnFdV1V+XUi6spWiY5hbOa23oOGB60LMBaIbJ8Nl0b68pfUaSf0h/sxzp68QkP6uq6thGFgkAjE1VVQuSnJNkbvp78y+S/Dj9ZzZVSd5YVdUjkpyW/s8BX0jynCRPSvJPSc7b+VRzk5xZVZWEDBpg2eIFaW9rHUyc7q5K/53Oly1eUGdZQBPp2QA0y2T4bDrmUHpn0Py89B/h/Mv0Hx19jyQHpT+MXp/+mhcluaCqqqc2vFoAYE9emOTgJDuS/F0p5SGllGVJHprk9zvHfDzJEUmeV0p5dinlv0sp3yqlfLqU8tQkr9k57n5JnlVz/TAtzWqpsnplR5IM+/A/8PPqlR2Z1TLapgEwDenZADTFZPhsujdHSr9w5/cbkvxFKeUbpZQ7Sim3llK+kmR5ktelv6Hun+R/qqp6TmPLBQD24IT070D+n1LK1wYmllJ+keR96f+M8fAkXymlnDvSE5RSPpDk0p0/2skMDbJiSXtOX7U0i9qGHsy4qK01p69amhVL2ptUGdAkejYATdPsz6Z7c03pY9PfMP+9lNJ195k77/j771VV/TjJ/yS5V5JPV1U1r5Ty0YZUCwDsScfO718aYd7/l+T9O//9hT08z7np7/0Pb0xZQNL/4f/4jkVZ39mVLVt7snBe/2mRjpCGGUnPBqCpmvnZdG9C6UN2fr9sd4NKKRdXVfXYJN/c+ZiPVFV1j1LK+3f3OACgIe658/v1I8y7YZd//3oPz3Plzu8Lx1sQMNSslirLjzyo2WUAzXfPnd/1bACaplmfTffm8h377fzeu6eBpZSrkjwmSWf6Tzl6b1VVa/a+PABgL92583vf3WeUUu4cYdxotu78vn8jigIAhtGzAZix9iaU/t+d3+8/lsGllM70B9M/T38w/Zaqqv5978oDAPbSQL8e7wXA7rHz+y3jfB4AYGR6NgAz1t5cvmNT+i/H8ZgkXxnLA0opN+68lMcFSY5J8sokD9vLGsdt246+nH3Ztbmu644ctuCAPHf54Zkze2/yeJgarOvMBNbzPbohyRFJ7jfK/H/a+X2kU4V3deTO779vRFEAwDB6NgAz1t6E0t9LcnySZ1ZV9ZqdNzbco1LKLVVVPSHJeUn+IskT977Mfffu8zbljEs607dLte867+c5+TGL88YTOkZ/IEwx1nVmAuv5mPw4/TuQl440s5Ty6TE+z3E7v1/TiKIAgGH0bABmrL05tOy8nd8PSfL0vVlIKeW29Afa/y/9l/KoxbvP25SPfXdoeJEkfSX52Hc78+7zNtVVCkwo6zozgfV8zH6U/l67zzuBq6qam+TEJCXJpQ2qCwAYSs8GYMYacyhdSrk8/UdL/ybJP+7tgnbeqOGpGeOlP8Zr246+nHFJ527HnHFJZ7btGHZPCZhSrOvMBNbzvXJe+ncEv7Cqqn3dEfzs9N806bdJvtGowgCAIfRsAGasvboIZynlsaWUB5ZSnrYvCyulbC+lnFhKaSmlzNqX5xirsy+7dtjRdHfXV/rHwVRmXWcmsJ6PXSllaynl2zu/xnSprRGe49OllMU7v37V6BoBAD0bgJlt2t4Z6rquOxo6DiYr6zozgfUcAAAApo9xhdJVVR2351HNcdiCAxo6DiYr6zozgfV8fCZzvwYA/kTPBmCmGO+R0hdUVfV3DamkwZ67/PC07OGqXC1V/ziYyqzrzATW83GbtP0aABhCzwZgRhhvKN2a5AtVVb28EcU00pzZLTn5MYt3O+bkxyzOnNnT9gomzBDWdWYC6/m4Tdp+DQAMoWcDMCOMd+v9B0lmJflgVVXvGeuDqqo6qqqqr4xz2Xv0xhM68uLHLh52dF1Llbz4sYvzxhM6JroEqIV1nZnAej4uk7pfAwCD9GwAZoRqH2/y2//gqmpN8tkkf5uk7Pz3P5VSdowyfnGStyd5dpKWUsqsMSxjfpLu7u7uzJ8/f5/q3LajL2dfdm2u67ojhy04IM9dfrij6ZiWrOvMBI1ez2+77ba0tbUlSVsp5baGFTqJTJV+DQCjmQn9Opn4nq1fAzDRxtqzxxVKJ0lVVVWSDyZ5efqb5oVJnl5K2brLmPYkb03ygiSzk1RJeksp+43h+TVNACbMDNrI1a8BmLJmSr9OJrZn69cATLSx9uxxH0JZ+v1zktftnPTEJN+tqmpRVVULqqp6X5JrkrwoyUCD/GKSo8e7bABgbPRrAJga9GwAZoLZjXqiUsq/V1X1uySfTvKwJJcnOTDJvPTvtS1J/ifJ20opGxu1XABg7PRrAJga9GwAprNGX2z2y0nOTn+DXJRk4HygryY5ppTyDM0SAJpOvwaAqUHPBmBaakgoXVXVrKqqTkryqyQvTP8e22rn7FuSrCml/KwRywIA9o1+DQBTg54NwHQ3rlB6Z6N8QZJfJvlYksPT3yh/muRlSW5OclCSi6uq+qvxlQoA7Av9GgCmBj0bgJlivEdK/yLJGUkWp79R/jLJs0spS0sppyc5Nsm16T/F6BtVVf3jOJcHAOw9/RoApgY9G4AZYbyh9BHpb5S/TXJSko5SyhcGZpZSfpXk0em/IcN+Sc6squpfxrlMAGDv6NcAMDXo2QDMCOMNpbckOTXJg0opZ5ZS+u4+oJTyv0kel2Rd+pvr26qq+nhVVY2+ySIAMDL9GgCmBj0bgBlh3EdKl1L+s5SybXeDSil3JFmZ5Mz0N80XJlk7zmUDAGOjXwPA1KBnAzAjjCuU3tkIxzq2t5RyUpK3p79prhjPsgGAsdGvAWBq0LMBmClqP72nlPK2JCcn6a172QDA2OjXADA16NkATEWzm7HQUsonq6q6sRnLBgDGRr8GgKlBzwZgqmnajRBKKec3a9kAwNjo1wAwNejZAEwl7s4LAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURigNAAAAAEBtZje7AAAAmEl6+0rWd3Zly9aeLJzXmmWLF2RWS9XssgCAXejXMLGE0gAAUJN1GzdnzdpN2dzdMzitva01q1d2ZMWS9iZWBgAM0K9h4rl8BwAA1GDdxs055ZwNQzZwk+Sm7p6ccs6GrNu4uUmVAQAD9Guoh1AaAAAmWG9fyZq1m1JGmDcwbc3aTentG2kEAFAH/RrqI5QGAIAJtr6za9gRV7sqSTZ392R9Z1d9RQEAQ+jXUB+hNAAATLAtW0ffwN2XcQBA4+nXUB+hNAAATLCF81obOg4AaDz9GuojlAYAgAm2bPGCtLe1phplfpWkva01yxYvqLMsAGAX+jXURygNAAATbFZLldUrO3Y7ZvXKjsxqGW0zGGDy6O0ruezXt+RrP7khl/36Fjd9Y9rYtV/fvSMP/KxfQ2PMbnYBAAAwE6xY0p4XPXZxzrikM7vmNy1VcvJjFmfFkvbmFQcwRus2bs6atZuG3Ayuva01q1d2eB9jWlixpD2nr1o6bD1fZD2HhhJKAwBADdZt3JyPf7czdz+esJTk49/tzDH3v5cNXWBSW7dxc045Z8Ow97Gbuntyyjkbcvqqpd7HmBZWLGnP8R2Lsr6zK1u29mThvP5LdjhCGhrH5TsAAGCC9faVrFm7aViQk2Rw2pq1m5wCD0xa3seYaWa1VFl+5EH5m4ffN8uPPEggDQ0mlAYAgAm2vrNryCnAd1eSbO7uyfrOrvqKAtgL3scAaCShNAAATLAtW0cPcvZlHEDdvI8B0EiuKQ0AABNs4bzWho4DqJv3MWaa3r7imtIwgYTSAAAwwZYtXpD2ttbdnvre3ta/wQswGQ28j93U3TPidaWrJIu8jzFNrNu4OWvWbhrSt9vbWrN6ZYebeUKDuHwHAABMsFktVZ529O43Yp92dLsjsIBJa1ZLldUrO5L0B9C7Gvh59coO72NMees2bs4p52wYtiP5pu6enHLOhqzbuLlJlcH0IpQGAIAJ1ttX8v/9dPcbsf/fTzent2+k4w8BJocVS9pz+qqlWdQ29BIdi9pac/qqpY4gZcrr7StZs3bTiGcDDExbs3aTfg0N4PIdAAAwwdZ3du320h1Jsrm7J+s7u7L8yINqqgpg761Y0p7jOxa51i7T0p76dYl+DY0ilAYAgAm2ZevuA+m9HQfQTLNaKoEc05J+DfVx+Q4AAJhgC+e17nnQXowDABpPv4b6CKUBAGCCLVu8IO1trcNuDjagStLe1n8KPADQHPo11EcoDQAAE2xWS5XVKzuSZNiG7sDPq1d2uCYrADSRfg31EUoDAEANVixpz+mrlmZR29BTfhe1teb0VUuzYkl7kyoDAAbo11APNzoEAICarFjSnuM7FmV9Z1e2bO3Jwnn9pwA74goAJg/9GiaeUBoAAGo0q6XK8iMPanYZAMBu6NcwsVy+AwAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNGx0CAECNevtK1nd2ZcvWniyc15plixdkVkvV7LIAxsz7GADjJZQGAICarNu4OWvWbsrm7p7Bae1trVm9siMrlrQ3sTKAsfE+BkAjuHwHAADUYN3GzTnlnA1Dgpwkuam7J6ecsyHrNm5uUmUAY+N9DIBGEUoDAMAE6+0rWbN2U8oI8wamrVm7Kb19I40AaD7vYwA0klAaAAAm2PrOrmFHFu6qJNnc3ZP1nV31FQWwF7yPAdBIQmkAAJhgW7aOHuTsyziAunkfA6CRhNIAADDBFs5rbeg4gLp5HwOgkYTSAAAwwZYtXpD2ttZUo8yvkrS3tWbZ4gV1lgUwZt7HAKan3r6Sy359S772kxty2a9vqe3eALNrWQoAAMxgs1qqrF7ZkZecs2HE+SXJ6pUdmdUyWtwD0FwD72OnnLMhVTLkhocD71zexwCmlnUbN2fN2k1D7hnQ3taa1Ss7smJJ+4Qu25HSAAAAwB6tWNKe01ctzaK2oZfoWNTWmtNXLZ3wAAOAxlm3cXNOOWfDsJvY3tTdk1PO2ZB1GzdP6PIdKQ0AABOst69kzdpNo86vkqxZuynHdyxylCEwqa1Y0p7jOxZlfWdXtmztycJ5/Zfs8N4FMHUMfDYd6UIdJfV8NhVKAwDABFvf2TXsKJRdlSSbu3uyvrMry488qL7CAPbBrJbKexXAFDYZPpu6fAcAAEywLVtH/9C/L+MAAGBfTYbPpkJpAACYYAvnte550F6MAwCAfTUZPpsKpQEAYIItW7wg7W2tGe2KfFX673S+bPGCOssCAGAGmgyfTYXSAAAwwWa1VFm9siNJhn34H/h59coONwoDAGDCTYbPpkJpAACowYol7Tl91dIsaht6GuSittacvmppVixpb1JlAADMNM3+bDp7Qp8dAAAYtGJJe47vWJT1nV3ZsrUnC+f1nxbpCGkAAOrWzM+mQmkAAKjRrJYqy488qNllAABA0z6bunwHAAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtZnd7AIAAAb09pWs7+zKlq09WTivNcsWL8islqrZZQEAu9CvARgvoTQAMCms27g5a9ZuyubunsFp7W2tWb2yIyuWtDexMgBggH4NQCO4fAcA0HTrNm7OKedsGLKBmyQ3dffklHM2ZN3GzU2qDAAYoF8D0ChCaQCgqXr7Stas3ZQywryBaWvWbkpv30gjAIA66NcANJJQGgBoqvWdXcOOuNpVSbK5uyfrO7vqKwoAGEK/BqCRhNIAQFNt2Tr6Bu6+jAMAGk+/BqCRhNIAQFMtnNfa0HEAQOPp1wA0klAaAGiqZYsXpL2tNdUo86sk7W2tWbZ4QZ1lAQC70K8BaCShNADQVLNaqqxe2ZEkwzZ0B35evbIjs1pG2wwGACaafg1AIwmlAYCmW7GkPaevWppFbUNP+V3U1prTVy3NiiXtTaoMABigXwPQKLObXQAAQNK/oXt8x6Ks7+zKlq09WTiv/xRgR1wBwOShXwPQCEJpAGDSmNVSZfmRBzW7DABgN/RrAMbL5TsAAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKjN7GYXAAAwYNuOvpx92bW5ruuOHLbggDx3+eGZM9s+dAAA6tXbV7K+sytbtvZk4bzWLFu8ILNaqmaXBQ3XrHVdKA0ATArvPm9TzrikM33lT9Pedd7Pc/JjFueNJ3Q0rzAAAGaUdRs3Z83aTdnc3TM4rb2tNatXdmTFkvYmVgaN1cx13aFHAEDTvfu8TfnYd4cG0knSV5KPfbcz7z5vU3MKAwBgRlm3cXNOOWfDkJAuSW7q7skp52zIuo2bm1QZNFaz13WhNADQVNt29OWMSzp3O+aMSzqzbUdfTRUBADAT9faVrFm7KWWEeQPT1qzdlN67H0kBU8xkWNeF0gBAU5192bXDjpC+u77SPw4AACbK+s6uYUeN7qok2dzdk/WdXfUVBRNgMqzrrikNADTVdV13NHQcADCx3ACO6WrL1tFDun0ZB5PVZFjXhdIAQFMdtuCAho4DACaOG8AxnS2c19rQcTBZTYZ13eU7AICmeu7yw1Pt4eCqquofBwA0T7NvigUTbdniBWlva81oH02r9O+EWbZ4QZ1lQcNNhnVdKA0ANNWslir77zdrt2MO2G+W04IBoIkmw02xYKLNaqmyemVHkgwL6wZ+Xr2yw+dSprzJsK4LpQGAplrf2ZU7tvXudswft/W6oQwANNFkuCkW1GHFkvacvmppFrUNvWzBorbWnL5qqcvUMG00e113TWkAoKkmw002AIDd06+ZSVYsac/xHYvc0JNpr5nrulAaAGiqyXCTDQBg9/RrZppZLVWWH3lQs8uACdesdd3lOwCAppoMN9kAAHZPvwagkYTSAEBTTYabbAAAu6dfA9BIQmkAoOmafZMNAGDP9GsAGsU1pQGAScENZQBg8tOvAWgEoTQAMGm4oQwATH76NQDj5fIdAAAAAADURigNAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURigNAAAAAEBthNIAAAAAANRmdrMLAAAY0NtXsr6zK1u29mThvNYsW7wgs1qqZpcFAAAwLTVrG0woDQBMCus2bs6atZuyubtncFp7W2tWr+zIiiXtTawMAABg+mnmNpjLdwAATbdu4+accs6GIR+GkuSm7p6ccs6GrNu4uUmVAQAATD/N3gYTSgMATdXbV7Jm7aaUEeYNTFuzdlN6+0YaAQAAwN6YDNtgQmkAoKnWd3YN2zu/q5Jkc3dP1nd21VcUAADANDUZtsGE0gBAU23ZOvqHoX0ZBwAAwOgmwzaYUBoAaKqF81obOg4AAIDRTYZtMKE0ANBUyxYvSHtba6pR5lfpvwP0ssUL6iwLAABgWpoM22BCaQCgqWa1VFm9siNJhn0oGvh59cqOzGoZ7SMTAAAAYzUZtsGE0gBA061Y0p7TVy3Norahp4ctamvN6auWZsWS9iZVBgAAMP00exts9oQ+OwDAGK1Y0p7jOxZlfWdXtmztycJ5/aeLOUIaAACg8Zq5DSaUBgAmjVktVZYfeVCzywAAAJgRmrUN5vIdAAAAAADURigNAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURigNAAAAAEBthNIAAAAAANRGKA0AAAAAQG2E0gAAAAAA1EYoDQAAAABAbYTSAAAAAADURigNAAAAAEBthNIAAAAAANRmdrMLqMO2HX05+7Jrc13XHTlswQF57vLDM2e2PJ7px7rOTGA9n956+0rWd3Zly9aeLJzXmmWLF2RWS9XssgAAmGF8LoWJNe1D6XeftylnXNKZvvKnae867+c5+TGL88YTOppXGDSYdZ2ZwHo+va3buDlr1m7K5u6ewWntba1ZvbIjK5a0N7EyAABmEp9LYeJN60PL3n3epnzsu0PDiyTpK8nHvtuZd5+3qTmFQYNZ15kJrOfT27qNm3PKORuGfPBPkpu6e3LKORuybuPmJlUGAMBM4nMp1GPahtLbdvTljEs6dzvmjEs6s21HX00VwcSwrjMTWM+nt96+kjVrN6WMMG9g2pq1m9J79z0SAADQQD6XQn2mbSh99mXXDjua7u76Sv84mMqs68wE1vPpbX1n17AjUXZVkmzu7sn6zq76igIAYMbxuRTqM21D6eu67mjoOJisrOvMBNbz6W3L1tE/+O/LOAAA2Bc+l0J9pm0ofdiCAxo6DiYr6zozgfV8els4r7Wh4wAAYF/4XAr1mbah9HOXH56WavdjWqr+cTCVWdeZCazn09uyxQvS3taa0V7iKv13O1+2eEGdZQEAMMP4XAr1mbah9JzZLTn5MYt3O+bkxyzOnNnT9k/ADGFdZyawnk9vs1qqrF7ZkSTDNgAGfl69siOz9rRnAgAAxsHnUqjPtN56f+MJHXnxYxcPO7qupUpe/NjFeeMJHc0pDBrMus5MYD2f3lYsac/pq5ZmUdvQUyEXtbXm9FVLs2JJe5MqAwBgJvG5FOpRlVKaXcNuVVU1P0l3d3d35s+fv0/PsW1HX86+7Npc13VHDltwQJ67/HBH0zEtWdeZCRq9nt92221pa2tLkrZSym0NK3SGaUS/TpLevpL1nV3ZsrUnC+f1nxrpSBQA9OvGaFS/hpnA51LYN2Pt2TMilAaA0djIbQz9GoCJpF83hn4NwEQba892CCUAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFAboTQAAAAAALURSgMAAAAAUBuhNAAAAAAAtRFKAwAAAABQG6E0AAAAAAC1EUoDAAAAAFCb2c0uYKxuu+22ZpcAwDSkvzSWvycAE0F/aSx/TwAmylh7TFVKmeBSxqeqqvsmub7ZdQAw7R1aSrmh2UVMVfo1ADXRr8dBvwagRrvt2VMhlK6SHJJka7NrAWDampfkxjLZm+Ikpl8DUAP9epz0awBqsseePelDaQAAAAAApg83OgQAAAAAoDZCaQAAAAAAaiOUBgAAAACgNkJpAAAAAABqI5QGAAAAAKA2QmmYAaqqWlhV1VOrqnp7VVXnV1V1c1VVZefXWc2uDwDQrwFgKqiq6pFVVb21qqoLqqq6vqqqu6qqur2qql9WVfWpqqqOa3aNMBXMbnYBQC1+3+wCAIA90q8BYBKrquq7SR4zwqw5SR648+v5VVV9JsnJpZRtddYHU4kjpWHm+W2SC5pdBACwW/o1AEw+h+z8fmOSDyZ5RpJlSZYneXWSG3bOf16Ss+ouDqYSR0rDzPD2JD9K8qNSyu+rqjo8SWdzSwIA7ka/BoDJ7eokb0ry5VJK793m/aCqqrOTXJrkQUn+T1VVHy2lfLfuImEqEErDDFBKWd3sGgCA3dOvAWByK6U8dQ/zb66q6jVJ1u6c9IwkQmkYgct3wCRWVdWBVVX9fucNjn5TVdV+o4zbv6qq7+8c11NV1WPrrhUAZir9GgAmvxr79UW7/PvIfa8YpjehNExipZQ/JvnXnT8uTvKPdx9TVVVLknPTfw2rviSrnB4EAPXRrwFg8quxX8/d5d93v8QHsJNQGia/j6b/ZkdJ8uYR9uaeluTvdv771aWUL9VVGAAwSL8GgMmvjn79uF3+/fN9eDzMCEJpmORKKXel/8ZHSXJ4kn8amLfzWlWv2Pnj+0spH6y3OgAg0a8BYCqY6H6980jrN+wy6Qv7VilMf0JpmBrOSvLLnf9+U1VV+1VV9awk79s57b+T/N9mFAYADDor+jUATHZnZeL69auSLNv57/8ppVy+z1XCNCeUhimglNKb5K07fzwsyX8l+UySKv03UXh+KaU0qTwAIPo1AEwFE9Wvq6p6XJL37PxxS5JTxl8tTF9CaZg6vpDkJzv/fVL6b55wZZK/K6Vsa1ZRAMAQ+jUATH4N7ddVVR2V5CtJZifpSfLMUsqWxpQK05NQGqaInXtqz9hl0k1JnlJK6W5SSQDA3ejXADD5NbJfV1W1OMkFSe6VpDfJs0sp321IoTCNCaVhiqiq6oFJ1uwy6cAkdzWpHABgBPo1AEx+jerXVVUdkuRbSQ5JUpK8oJTytYYUCdOcUBqmgKqqFiZZl+TgJLfsnDwvyRubVhQAMIR+DQCTX6P6dVVVByf5f0mO2DnpFaWUzzSqTpjuhNIwyVVVdWCSb6S/0d2e5PgkX905+6VVVd23SaUBADvp1wAw+TWqX1dV1Zbkm0k6dk56QynlI42tFqY3oTRMYlVVzU7yxSSPTLIjyTNKKVckeVv6Tw1qzZ/uGgwANIF+DQCTX6P6dVVVB6Q/2F66c9K7SinvnYiaYToTSsPk9tEkT9n57xeXUr6ZJKWUn6b/zr5J8oKqqo5sRnEAQBL9GgCmgnH366qq5uwce+zOSR8spbxlguqFaa3qv+EoMNlUVfW2JKt3/rimlPK2u81/aJKfJqmSfLaU8pzdPNdxSR6wy6SDk7xv578vTfKJXceXUs4aR+kAMGPo1wAw+TWqX1dV9eUkT9/544VJXpn+o6xHs62U8st9LhymMaE0TEJVVb0wf9rwPLOU8sJRxn0hyTOT9CV5eCnlylHGnZXkH8e6/FJKtVcFA8AMpF8DwOTXyH5dVdXehmjXlVIO38vHwIzg8h0wyVRVdUL6TytK+m+c8OLdDF+T/obZkuSdE1waALCTfg0Ak59+DZOXI6UBAAAAAKiNI6UBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACojVAaAAAAAIDaCKUBAAAAAKiNUBoAAAAAgNoIpQEAAAAAqI1QGgAAAACA2gilAQAAAACozf8PKL3NDGCv2bEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1800x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 为在Matplotlib中显示中文，设置特殊字体\n",
    "plt.rcParams[\"font.sans-serif\"] = [\"SimHei\"]\n",
    "# 正确显示负号\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "plt.rcParams.update({'font.size': 20})\n",
    "# 创建一个图形框\n",
    "fig = plt.figure(figsize=(18, 6), dpi=100)\n",
    "# 计算相关系数\n",
    "corr = data.corr(method='pearson')\n",
    "var = [('x1', 'x2'), ('x1', 'x3'), ('x2', 'x3')]\n",
    "for i, (v1, v2) in enumerate(var):\n",
    "    ax = fig.add_subplot(1, 3, i + 1)\n",
    "    ax.set_xlabel(f'${v1}$')\n",
    "    ax.set_ylabel(f'${v2}$')\n",
    "    ax.scatter(data[v1], data[v2])\n",
    "    ax.annotate(f'相关系数：{corr[v1][v2]:.2f}', (0.5, 0.9),\n",
    "                xycoords='axes fraction', ha='center', va='center')\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "plt.savefig('multicollinearity_data.png', dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(X, Y):\n",
    "    '''\n",
    "    训练线性回归模型\n",
    "    '''\n",
    "    model = sm.OLS(Y, X)\n",
    "    re = model.fit()\n",
    "    return re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.431\n",
      "Model:                            OLS   Adj. R-squared:                  0.396\n",
      "Method:                 Least Squares   F-statistic:                     12.14\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00307\n",
      "Time:                        17:46:01   Log-Likelihood:                -23.867\n",
      "No. Observations:                  18   AIC:                             51.73\n",
      "Df Residuals:                      16   BIC:                             53.51\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.9094      0.360     -2.525      0.023      -1.673      -0.146\n",
      "x1             0.9719      0.279      3.484      0.003       0.380       1.563\n",
      "==============================================================================\n",
      "Omnibus:                        6.365   Durbin-Watson:                   2.636\n",
      "Prob(Omnibus):                  0.041   Jarque-Bera (JB):                1.638\n",
      "Skew:                           0.020   Prob(JB):                        0.441\n",
      "Kurtosis:                       1.523   Cond. No.                         2.92\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "# 用x1搭建回归模型\n",
    "Y = data['y']\n",
    "X = sm.add_constant(data['x1'])\n",
    "re = train_model(X, Y)\n",
    "print(re.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.566\n",
      "Model:                            OLS   Adj. R-squared:                  0.539\n",
      "Method:                 Least Squares   F-statistic:                     20.84\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):           0.000318\n",
      "Time:                        17:46:01   Log-Likelihood:                -21.442\n",
      "No. Observations:                  18   AIC:                             46.88\n",
      "Df Residuals:                      16   BIC:                             48.66\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.1755      0.315      3.734      0.002       0.508       1.843\n",
      "x2            -1.1130      0.244     -4.565      0.000      -1.630      -0.596\n",
      "==============================================================================\n",
      "Omnibus:                        6.451   Durbin-Watson:                   0.642\n",
      "Prob(Omnibus):                  0.040   Jarque-Bera (JB):                1.649\n",
      "Skew:                           0.035   Prob(JB):                        0.438\n",
      "Kurtosis:                       1.519   Cond. No.                         2.92\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "# 用x2搭建回归模型\n",
    "X1 = sm.add_constant(data['x2'])\n",
    "re1 = train_model(X1, Y)\n",
    "print(re1.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.997\n",
      "Model:                            OLS   Adj. R-squared:                  0.997\n",
      "Method:                 Least Squares   F-statistic:                     2488.\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):           1.21e-19\n",
      "Time:                        17:46:01   Log-Likelihood:                 23.319\n",
      "No. Observations:                  18   AIC:                            -40.64\n",
      "Df Residuals:                      15   BIC:                            -37.97\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.2036      0.034      5.952      0.000       0.131       0.277\n",
      "x1             0.9719      0.021     46.398      0.000       0.927       1.017\n",
      "x2            -1.1130      0.021    -53.133      0.000      -1.158      -1.068\n",
      "==============================================================================\n",
      "Omnibus:                        0.411   Durbin-Watson:                   2.369\n",
      "Prob(Omnibus):                  0.814   Jarque-Bera (JB):                0.514\n",
      "Skew:                           0.278   Prob(JB):                        0.774\n",
      "Kurtosis:                       2.387   Cond. No.                         4.26\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "# 用不相关的x1，x2搭建回归模型\n",
    "X2 = sm.add_constant(data[['x1', 'x2']])\n",
    "re2 = train_model(X2, Y)\n",
    "print(re2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.431\n",
      "Model:                            OLS   Adj. R-squared:                  0.396\n",
      "Method:                 Least Squares   F-statistic:                     12.14\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00307\n",
      "Time:                        17:46:01   Log-Likelihood:                -23.867\n",
      "No. Observations:                  18   AIC:                             51.73\n",
      "Df Residuals:                      16   BIC:                             53.51\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.9094      0.360     -2.525      0.023      -1.673      -0.146\n",
      "x1             0.9719      0.279      3.484      0.003       0.380       1.563\n",
      "==============================================================================\n",
      "Omnibus:                        6.365   Durbin-Watson:                   2.636\n",
      "Prob(Omnibus):                  0.041   Jarque-Bera (JB):                1.638\n",
      "Skew:                           0.020   Prob(JB):                        0.441\n",
      "Kurtosis:                       1.523   Cond. No.                         2.92\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# 用强相关的x1，x3搭建模型\n",
    "Y = data['y']\n",
    "X = sm.add_constant(data['x1'])\n",
    "re = train_model(X, Y)\n",
    "print(re.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.483\n",
      "Model:                            OLS   Adj. R-squared:                  0.451\n",
      "Method:                 Least Squares   F-statistic:                     14.96\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00136\n",
      "Time:                        17:46:01   Log-Likelihood:                -23.007\n",
      "No. Observations:                  18   AIC:                             50.01\n",
      "Df Residuals:                      16   BIC:                             51.79\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.4717      0.452     -3.254      0.005      -2.431      -0.513\n",
      "x3             1.0901      0.282      3.867      0.001       0.493       1.688\n",
      "==============================================================================\n",
      "Omnibus:                        4.131   Durbin-Watson:                   2.709\n",
      "Prob(Omnibus):                  0.127   Jarque-Bera (JB):                1.375\n",
      "Skew:                           0.060   Prob(JB):                        0.503\n",
      "Kurtosis:                       1.652   Cond. No.                         4.41\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "X1 = sm.add_constant(data['x3'])\n",
    "re1 = train_model(X1, Y)\n",
    "print(re1.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.484\n",
      "Model:                            OLS   Adj. R-squared:                  0.415\n",
      "Method:                 Least Squares   F-statistic:                     7.040\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00698\n",
      "Time:                        17:46:01   Log-Likelihood:                -22.989\n",
      "No. Observations:                  18   AIC:                             51.98\n",
      "Df Residuals:                      15   BIC:                             54.65\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.5437      0.622     -2.480      0.025      -2.870      -0.217\n",
      "x1            -0.1677      0.959     -0.175      0.864      -2.213       1.877\n",
      "x3             1.2604      1.017      1.240      0.234      -0.907       3.428\n",
      "==============================================================================\n",
      "Omnibus:                        3.209   Durbin-Watson:                   2.726\n",
      "Prob(Omnibus):                  0.201   Jarque-Bera (JB):                1.239\n",
      "Skew:                           0.079   Prob(JB):                        0.538\n",
      "Kurtosis:                       1.724   Cond. No.                         14.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X2 = sm.add_constant(data[['x1', 'x3']])\n",
    "re2 = train_model(X2, Y)\n",
    "print(re2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "检测假设x1和x3同时不显著：\n",
      "<F test: F=7.040116180053257, p=0.006977508919044013, df_denom=15, df_num=2>\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>VIF Factor</th>\n",
       "      <th>features</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>11.049590</td>\n",
       "      <td>const</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>12.907932</td>\n",
       "      <td>x1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.061102</td>\n",
       "      <td>x2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>12.969034</td>\n",
       "      <td>x3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   VIF Factor features\n",
       "0   11.049590    const\n",
       "1   12.907932       x1\n",
       "2    1.061102       x2\n",
       "3   12.969034       x3"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 检测多重共线性\n",
    "print('检测假设x1和x3同时不显著：')\n",
    "print(re2.f_test(['x1=0', 'x3=0']))\n",
    "vif = pd.DataFrame()\n",
    "X = sm.add_constant(data[['x1', 'x2', 'x3']])\n",
    "vif['VIF Factor'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "vif['features'] = X.columns\n",
    "vif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+wAAAIMCAYAAABv6ia+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVe0lEQVR4nO3de5xcdX34/9c7AVmBZDGEdFehmOJ1DWIjRlOlghYb0MjPFrXV0NYbNtRLq/WCl2+MN6q2tRY1VVpLJbVqL18lFiK1XBQB88UoGKMPAQMCbgCJbiKwCLvv3x9nJkw2s5fZnd05M/N6Ph7nMcxnPufMew6bOec953Pen8hMJEmSJElSucxrdQCSJEmSJGl/JuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJXRAqwOYLRERwCOBPa2ORZKkigXATzMzWx1IJ/BYL0kqqaYd7zs2Yac4gN/W6iAkSRrjSOD2VgfRITzWS5LKqinH+05O2PcA3HrrrSxcuLDVsUiSutzu3bs56qijwKvBzeSxXpJUKs0+3ndywg7AwoULPYhLktTBPNZLkjqVReckSZIkSSohE3ZJkiRJkkrIhF2SJEmSpBIyYZckSXtFxLyIeHdE3BoRD0TEbRGxdhrbOT0itkfEvRFxYUQcUafPiRFxbaXP1yPimOZ8CkmSOoMJuyRJqvWeyvI14CzgBuCTEfHiqW4gIk4DvgjcDZwN/HrleW2fpwGbgYOAdwAJXBQRB834E0iS1CFKmbBHxMERcV5EDEbEPRGxNSJe0uq4JEnqZJWr4G8F3pSZr8jM84BTgDuAV01xG/OAjwLbgJMz82PAqcDKiFhd0/UjwC+AZ2fm31X6LATObM6nkSSp/ZUyYQc+BryQ4oD/ZuBnwBfGHOglSVJzHQKsBz5RbcjMYeBGYL8h7eM4FlgKfKKyLpn5U+Bi4DSAiHgEcALwz5m5q9LnHuDz1T6SJKmECXtEHAj8MfC2zPxwZv4DsAr4MfDylgYnSVIHy8ybM/OczHyw2lY5Lh8LfH+Km3lS5fHKMe3XA8dV/vuJFOcgE/XZT0QcFBELqwuwYIoxSZLUlkqXsAOLgAOBe6sNmTkKPAjc36qgJEnqUq+mGKr+2Sn2P6zyePOY9ruAo6bQZ3FE9Iyz7bOBoZrltinGJElSWypdwp6Zd1D8iv/OiHh8RCyIiLcBjwP+b2ujkySpe0TEYooCdJdl5temulrl8b4x7cM8dEW82ufeOn1g/Cvn5wC9NcuRU4xJkqS2dECrAxjHCyiGyf2w8jyBN2bml8ZboVJVtrayrMPkJEmamQ0U97U3UghumOK4PQ8YrWkPHjpOD9e0MaYP7Hs83ysz76dmtF3E2NUlSeospbvCXvHnQD/w38AFFFVkz46IFROs4zA5SZKaJCLOBE4HXpeZNzaw6k6KxPtRY9oXA3tq+sBDQ+Rr+1DTT5Kkrla6hD0ingK8EXhZZr4gM/8IeAJwD/DPE6zqMDlJkpogIp4O/D3wmcw8v8HVv1d5XDmmfTkPJeo3UgyZr9fnvswcavA9JUnqSKVL2IGTgF2Z+YVqQ2beSVHsZqAyFcx+MvP+zNxdXfDXeUnSDIyMJlffdDdf/u7tXH3T3YyMZqtDmhMR8RhgE8U86n/W6PqZ+RPgOuC1URmzHhGHU8znfnmlz/3AJcArqgXmKtXoXwpcMfNPIUnS5NrhWF/Ge9gDOCAiIjNr99jBlccyxixJ6iCbtw2yftN2BoeG97b19/awbvUAq5b1tzCyOXEBxZzrHwVOr71PPDM3RsRK4FGZ+R8TbOMcijnVPx0R/wn8H6AHOK+mz4eBbwBfiIgNwFkUQ+Rf18TPIklSXe1yrI99c+LWi4jfB/4DWJ2ZX6m0PRK4lqJ4zVE5haAr87MODQ0NsXDhwtkMWZLUQTZvG2Ttxq2MPdBU09YNa5ZP60C+e/duent7AXorI8FKp3Il/GfjvZ6ZERHnAy/IzMXj9ats663AByh+aL8feENmfnpMnzXAP1AUthsF3p+Z6xqI12O9JKlhs3Wsh+Yf78uYsB8C/ICiWM01FAf5p1NcYT8zM8+bYPXa7XgQlyQ1ZGQ0edaHLt3n1/ZaAfT19nDl257D/HmNVShvh4S92SLiUcBxwHWZefs4fQ4HVgA3NFjczmO9JKlhs3msh+Yf70t3D3tm3gOcSDHn+hOAE4A7KarA/2PrIpMkdbotO3aNewCHYq6ywaFhtuzYNXdBtbHMvD0zLxovWa/0uTszL240WZckaTra7VhfyvvBM/PHFFPJSJI0Z+7cM/4BfDr9JElSubTbsb50V9glSWqVJQt6mtpPkiSVS7sd603YJUmqWLF0Ef29PYx3x1pQVJBdsXTRXIYlSZKapN2O9SbskiRVzJ8XrFs9ALDfgbz6fN3qgWkVoZEkSa3Xbsd6E3ZJkmqsWtbPhjXL6evddyhcX2/PjKZ5kSRJ5dBOx/rSTevWLE71IkmaiZHRZMuOXdy5Z5glC4qhcTP5tb0bp3WbbR7rJUkz0exjPTT/eF/KKvGSJLXa/HnBymMOb3UYkiRplrTDsd4h8ZIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZL2ExFHRMStEXFig+u9JyJyvKWm38vG6bOm2Z9FkqR2dUCrA5AkSeUSEQ8HvggcOY3V/wu4sU77a4FH1jx/KvAD4INj+n1zGu8pSVJHMmGXJEl7RcRi4ELg6Omsn5nXA9eP2eYS4FPAWTXNxwOXZubGaYYqSVLHK+WQ+Ig4caLhdI0Oz5MkSVP2cuBB4IVN3ObbgJ8AGwEiYh6wHLiqie8hSVLHKesV9h8AZ9Rp/23glcAtcxuOJEld40Lg48BRzdhYRCwC/hRYm5kjlebHA4cCr4mIc4EDgcuAt2Tmj5rxvpIkdYJSJuyZeQeVX+FrRcQbgQsyc8fcRyVJUuerHmMjolmbfBVwL/CFmranVx4DWAccBrwBuCIinpiZv6i3oYg4CDiopmlBs4KUJKmMSpmw1xMRq4HjgBe3OhZJkjS5ytD3tcA/Zub9NS/9iCJB/3hmZqXvRcC1wBqKK/z1nE2R4EuS1BVKeQ/7ON4JbMzMm1sdiCRJmpLnAkuBz9Q2ZuZVmXluNVmvtG0Fvg+cNMH2zgF6a5bpVLGXJKlttMUV9og4nmL43FkT9HGYnCRJ5fJS4HuZecMU+98L9I33YuUq/d4r9U0cti9JUim1yxX21wPXVH59H8/ZwFDNcttcBCZJkvYXEQcCL6KYl33sax+JiDVj2g4AngD8bG4ilCSp/EqfsEfEAuAljBlOV4fD5CRJKo9nAouAzXVeWwq8vTI6ruotwELgq3MQmyRJbaEdhsSfRjHdy5cm6uQwOUmSZl9EnAwckJkXT9L1OcAwUG903F8DVwJXRcQVwOOA5wPfZvIf6CVJ6hqlv8JOcf/blZl5V6sDkSRJvBN43xT6nQR8OzN/NfaFzLwGOJViWrezgGUUSfxJmTncxFglSWprpb7CHhE9wO8A61sdiyRJ3aQyK8t+w9Uy88Qprn/CJK9fAlwyndgkSeoWZb/C/ltAD3BVqwORJEmSJGkulT1hPwl4APh/rQ5EkiRJkqS5VOqEPTPfnZkPy8z7Wh2LJEmSJElzqdQJuyRJkiRJ3cqEXZIkSZKkEjJhlyRJkiSphEzYJUmSJEkqIRN2SZIkSZJKyIRdkiRJkqQSMmGXJEmSJKmETNglSZIkSSohE3ZJkiRJkkrIhF2SJEmSpBIyYZckSZIkqYRM2CVJkiRJKiETdkmSJEmSSsiEXZIkSZKkEjJhlyRJkiSphEzYJUmSJEkqIRN2SZIkSZJKyIRdkiRJkqQSMmGXJEmSJKmETNglSZIkSSohE3ZJkiRJkkrIhF2SJEmSpBIyYZckSZIkqYRM2CVJkiRJKiETdkmSJEmSSsiEXZIkSZKkEjJhlyRJkiSphEzYJUmSJEkqIRN2SZIkSZJK6IBWByBJ6kwjo8mWHbu4c88wSxb0sGLpIubPi1aHJUmS1DZM2CVJTbd52yDrN21ncGh4b1t/bw/rVg+wall/CyOTJElqHw6JlyQ11eZtg6zduHWfZB1g59AwazduZfO2wRZFJkmS1F5M2CVJTTMymqzftJ2s81q1bf2m7YyM1ushSZKkWibskqSm2bJj135X1mslMDg0zJYdu+YuKEmSpDZlwi5Japo794yfrE+nnyRJUjczYZckNc2SBT1N7SdJktTNTNglSU0xMpqMZnLYww8ct09QVItfsXTR3AUmSZLUpkzYJUkztnnbIM/60KW8/B+/xS/ue6Bun+oM7OtWDzgfexuIiCMi4taIOHEa674sIrLOsmZMv9MjYntE3BsRF0bEEc2KX5KkTuA87JKkGalO4zZZ3fc+52FvGxHxcOCLwJHT3MRTgR8AHxzT/s2a9zit8h7fBD4FvKLy/KRpvqckSR2n9Al7RKwCLgb+IDO/0Op4JEkPmWgat6rDDj6QT/zhcp5xzOFeWW8DEbEYuBA4egabOR64NDM3jvMe84CPAtuAkzNzOCL+HfhxRKzOzE0zeG9JkjpGqYfER8TBwAbgf03WJal8JpvGDeAX9z7AvHlhst4+Xg48CLxwOitXkvHlwFUTdDsWWAp8IjOHATLzpxQ/0J82nfeVJKkTlTphB94DPAp4XYvjkCTV4TRuHelCimHpd09z/ccDhwKviYi7I2J3RHw5Ih5X0+dJlccrx6x7PXDceBuOiIMiYmF1ARZMM0ZJktpCaRP2iHgy8BfARzPzh62OR5K0P6dx6zyZuSMzR2awiadXHgNYB3wYWAlcERGHVV6rPt48Zt27gKMm2PbZwFDNctsM4pQkqfRKmbBHRACfBn4JXBkRp0bE4S0OS5I0xoqli+jv7WG8we5O49aVfgS8ATgpMz+eme8HVgG/BlSrxFf/ZO4bs+4wE181PwforVmmWxRPkqS2UMqEHfgjil/oDwM+A/wncEtEvHy8FRwmJ0lzb/68YN3qAYD9knancetOmXlVZp6bmVnTthX4Pg9VgB8Gkv3PQwI4aIJt35+Zu6sLsKe50UuSVC6lS9grV9f/D/AA8PzMPAJYAnwZ+JeIWD7Oqg6Tk6QWWLWsnw1rltPXu++w977eHjasWe40bqq6F+ir/PdOiuT8UWP6LMYkXJKkvco4rdvjgd8A/ikzLwLIzD0RcRbwh5Vla531zgH+tub5AkzaJWlOrFrWz8kDfWzZsYs79wyzZEExDN4r690nIj4CXFc7pVtEHAA8Abi80vS9yuNK4Jaa1ZdTJPOSJIlyJuzVGx2vqG3MzKGIuIv9f42vvn4/cH/1eXGhXpI0V+bPC1YeY7kRsRQ4JSL+vXJsBngLsBD4KkBm/iQirgNeGxFfyMys1Ko5BbigJVFLklRCZUzYByuPWdtY+XX+cIoKspIkqQUi4mTggMy8eJwuf00xXdtVEXEF8Djg+cC3KerSVJ0DfB74dET8J8XtcD3AebMVuyRJ7aZ097Bn5g6K4XEvHvPSHwHzgUvnPChJklT1TuB9472YmdcAp1Lco34WsIwiiT8pM4dr+n0BeBvwJ8DFFMPhz6oUqJMkSUDUFHEtjUo1+AuArwAXUdz3dhbwHWBlZo5OYRsLgaGhoSEWLlw4m+FKkjSp3bt309vbC9BbqXAuICIeBRxHcd/77Q2u67FeklQqzT7el3FIPJn5rxGxm+Ketw9RVIz/d+AvppKsS5Kk9lBJ0htK1CVJ6halTNgBMnMTsKnVcUiSJEmS1Aqlu4ddkiRJkiSZsEuSJEmSVEom7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJXRAqwOQJElqpZHRZMuOXdy5Z5glC3pYsXQR8+dFq8OSJMmEXZIkda/N2wZZv2k7g0PDe9v6e3tYt3qAVcv6WxiZJEkOiZckSV1q87ZB1m7cuk+yDrBzaJi1G7eyedtgiyKTJKlgwi5JkrrOyGiyftN2ss5r1bb1m7YzMlqvhyRJc8OEXZIkdZ0tO3btd2W9VgKDQ8Ns2bFr7oKSJGkME3ZJktR17twzfrI+nX6SJM0Gi85JUgex2rU0NUsW9DS1nyRJs8GEXZI6hNWupalbsXQR/b097BwarnsfewB9vcWPXpIktYpD4iWpA1jtWs0WEUdExK0RceI01j08Iv4lInZFxP0RcX1E/PaYPi+LiKyzrGnWZ5jI/HnButUDRSxj4688rls94AgVSVJLmbBLUpuz2rWaLSIeDnwROHIa684D/gN4EfBJ4E3AAuDiiDi6putTgR8AZ4xZvjmj4Buwalk/G9Ysp69332Hvfb09bFiz3JEpkqSWc0i8JLW5Rqpdrzzm8LkLTG0pIhYDFwJHT9Z3HC8AngWszMxrK9v8JvAd4A+Bv6r0Ox64NDM3zizimVm1rJ+TB/qs/SBJKqVSJuwR8Ujg9jov/VNmvnqu45GkMrPatZrs5cCDwAuBa6ex/t3Aa6vJesW2yuMRsPcq/HLgUzOIs2nmzwt/zJIklVIpE3aKYXIjwCuB0Zr2G1sTjiSVl9Wu1WQXAh8HjprOypn5TfYf1v7UyuP3K4+PBw4FXhMR5wIHApcBb8nMH03nfSVJ6kRlTdiPB36YmZ9tdSCSVHZWu1YzZeYOgIimDgl/E7Ab+FLl+dMrjwGsAw4D3gBcERFPzMxf1NtIRBwEHFTTtKCZQUqSVDZlLTp3PHBVq4OQpHZgtWuVWUQ8G3gx8DeZuavS/COKBP2kzPx4Zr4fWAX8GjBRlfizgaGa5bZZC1ySpBIoc8K+PCJuiYj7IuKaiDhlohUi4qCIWFhd8Fd3SV3Eatcqo4g4FPhn4Hs8VGyOzLwqM8/NzKxp20oxZP6kCTZ5DtBbszRcxV6SpHZSuiHxEfFoYAmwB9gA3AO8GvhKRJyQmeNdeT+bYlidJHUlq12rhD5FcdX8BZn5qyn0vxfoG+/FzLwfuL/6vMnD9iVJKp3SJewUxebeCWzIzJ8DRMT5wHaKe+DGS9jPAf625vkCHConqctY7VplERF/DrwM+KPM3D7mtY8A19VO6RYRBwBPAC6fwzAlSSq10iXsmXkr8MExbXsi4ksUB/7x1vNXd0mSSiAiXgT8DfCJzLygTpelwCkR8e+V4zfAW4CFwFfnKExJkkqvdAn7BO4FFkXEgZn5QKuDkSSpG0XEycABmXnxOK8vobhvfQj4dkTUFpG7IzP/B/hr4Ergqoi4Angc8Hzg28BnZjN+SZLaSekS9og4A1iWmW8b89KxwG6TdUmSWuqdFHOo103YgWdRFISD/ZPvK4D/ycxrIuJUikJ0ZwE7KZL492bmcPNDliSpPZUuYac4yL8hIj6dmTcBRMSzgN8F/rOlkUmS1CUy82b2nymQzDxxkvX+q956dfpdAlwyzfAkSeoKZZzW7bPAz4GvR8S5lYJzX6UYWveuVgYmSZIkSdJcKV3Cnpm7gWcD1wGvpLin7UJgRWb+qJWxSZIkSZI0V8o4JJ7MvAE4tdVxSJIkSZLUKqW7wi5JkiRJkkp6hV2SJEmSVG4jo8mWHbu4c88wSxb0sGLpIubPm7TuqBpgwi5JkiRJasjmbYOs37SdwaGHZuPs7+1h3eoBVi3rb2FkncUh8ZIkSZKkKdu8bZC1G7fuk6wD7BwaZu3GrWzeNtiiyDqPCbskSSUVEadHxN9FxCsj4oAxr/13q+KSJHWvkdFk/abtZJ3Xqm3rN21nZLReDzXKhF2SpBKKiNcBHwcOBt4CfDMiFtV0OaElgUmSutqWHbv2u7JeK4HBoWG27Ng1d0F1MBN2SZLK6XXA72bmmcCxwHeBS2uSdqv6SJLm3J17xk/Wp9NPEzNhlySpnPoz8zqAzHwwM18LXApcFhGHQ93RiJIkzaolC3qa2k8TM2GXJKmcfhYRS2sbMvNNwGWVxZleJElzbsXSRfT39ow7zCsoqsWvWLponB5qhAm7JJXAyGhy9U138+Xv3s7VN91toRYB/C/wJ2MbM/PPgcsBL11Ikubc/HnButUDwP73ZlWfr1s94HzsTRKZnXlSGBELgaGhoSEWLlzY6nAkaVzOY9oddu/eTW9vL0BvZu6erH9EPAw4IDPvHef1X8/MnzQ5zLbisV6SWsfzl/oaPd5PxoRdklqoOo/p2G/i6m/SG9Ys7+qDXidp9gFcHuslqRVGRpMtO3Zx555hFh96ECT87J77WbKgGAbf7VfWm3289/43SWqRyeYxDYp5TE8e6Ov6g5/2FRGPA54ELKH4c7kL2JaZN7Q0MElSR5voqvrKYw5vYWSdy4RdklqkkXlMPQgqIp4I/ClwOtBXba48ZqXPHcAXgU9l5g/mPEhJUscab1TgzqFh1m7c6qjAWWLCLkkt4jymmoqIOAb4EPAi4D7gG8CngJuAuymS9kXAY4BnAK8GXh8R/wW8LTN/3Iq4JUmdw1GBrWPCLkkt4jymmqLtwPcoKsb/V2beM1HniDiE4ir8Gyvr+gckSZoRRwW2jtO6SVKLOI+pJhIRI5X/fHFmHp+ZF0yWrANk5j2Z+S+ZuRx46exGKUnqBo4KbB0TdklqEecx1SQCIDMvnO4GMvPLzQtHktStHBXYOibsktRCq5b1s2HNcvp69z3A9fX2WLxFnTnvqiSp7TgqsHW8h12SWmzVsn5OHujbO6ep85hKkqQyqY4KXLtxK8G+vyg7KnB2eYVdkkpg/rxg5TGHc9pTHsXKYw73gKf9RMTiiPjLiPhARDyzpv1tEXFLROyJiEsi4thWxilJ6kyOCmwNr7BLklRyEdEHXAs8stL09ohYAywG1gPfAX4OPAe4PCKOy8zbWhKsJKljOSpw7pmwS5JUfm8BeoGXAD8B/gF4P/AL4GmZ+T2AiDgJ2Az8JfDnrQhUktTZqqMCNTccEi9JUvmdCnwmM/8jM7cAbweWAhdWk3WAzLwM+DzwvNaEKUmSmsmEXZKk8vt14Hs1z79fefxOnb7fBo6e9YgkSdKsM2GXJKn8HmDf29iGK4+/rNN3GMadeUeSJLURE3ZJksrvp8CRNc9/Cbwe+GGdvkcDd81FUJIkaXZZdE6SpPLbCjyj+iQz7wc+MU7f5wDXzUVQkiRpdpmwS5JUfu9l3yvsdUXErwGDwOdmPSJJkjTrTNglqYlGRtO5SdV0mfkj4EcAEfGSzPziOP3uAH5vLmOTJEmzx4Rdkppk87ZB1m/azuDQ8N62/t4e1q0eYNWy/hZGpg7zuYh4RGZ+qtWBSJKk2WXROUlqgs3bBlm7ces+yTrAzqFh1m7cyuZtgy2KTB3ofOCTEfHOei9GxMqI+MbchiRJkmaDCbskzdDIaLJ+03ayzmvVtvWbtjMyWq+H1JjMfDXwYeB9EfHRantEPD4i/gu4EvjNVsUnSZKap6Eh8RFxOvAs4Hrgs5n5YM1r/52Zz29yfJJUelt27NrvynqtBAaHhtmyYxcrjzl87gJTx8rMsyPiDuBvIuIIimneXknx5/Yp4H2tjE+SJDXHlK+wR8TrgI8DBwNvAb4ZEYtqupzQ5NgkqS3cuWf8ZH06/aQpOg+4DHgZ8Brg34EnZuZZmek9GJIkdYBGhsS/DvjdzDwTOBb4LnBpTdJuGWRJXWnJgp6m9pMmEhEHRsQbgZuAk4DvUFxZPxC4tZWxSZKk5mokYe/PzOsAMvPBzHwtcClwWUQcDnVv35Skjrdi6SL6e3vG/dUyKKrFr1i6aJweUkNuAP4WuAt4YWY+FXg58ELgoog4tBlvEhFHRMStEXHiNNc/KyJ+HBH3RMS/RMTBdfqcHhHbI+LeiLiwMrxfkiRVNJKw/ywiltY2ZOabKIbjXYZTxEnqUvPnBetWDwD7DzWqPl+3esD52NUs8ymGwB+Xmf8NkJlfoEjYn0HxQ/qMEt+IeDjwReDIaa7/euATwPeBdwG/DXxyTJ/TKu9xN3A28OuV55IkqaKRhP1/gT8Z25iZfw5cDjjWU1LXWrWsnw1rltPXu+9XYV9vDxvWLHcedjXTYzPzM5k5WtuYmZcAvwP8BkWl+GmJiMUUx/zHTXP9hcD7ga8Cp2XmR4HTgT+KiCdX+swDPgpsA07OzI8BpwIrI2L1dGOXJKnTNHJV/HXj9c/MN0TEXzcnpP1FxGuAT2eml6ckldaqZf2cPNDHlh27uHPPMEsWFMPgvbKuaar7h5OZ41YvzMxvRcQJwOYZvO/LgQcprthfO431TwIWAh+t/qiQmd+OiOsq27yeohbOUuBPq58nM38aERcDpwGbZhC/JEkdY8oJe2b+CvjVBK//pCkRjRERfRTzzUpS6c2fF07dpqbIzEZGwdWutz0injmDt76QYlaYo6a5/pMo6tp8c0z79cBxNX1g/5EA11Ncaa8rIg4CDqppWjDNGCVJagvTOhmAYshcRPxlRHyg9sQgIt4WEbdExJ6IuCQijp1hjOcCvTPchiRJXSMzp10tPjN3ZObIDN7+MODuzPzlmPa7eOhHgMMqjzdP0Kees4GhmuW2GcQpSVLpTatQXOWq97XAIytNb4+INcBiYD3FFDM/B54DXB4Rx2VmwwfVyn1svw98BnjVdGKVJElzKoB767QP89AV8epw//sm6FPPORQV8qsWYNIuSepg073C/haKq94voahIex1FgZk/AZ6WmSsz8ynAycChwF82+gYRsYCiouwGplA8JyIOioiF1QWHyUmS1ArD1L//PnhoOPswxbD5sechtX32k5n3Z+bu6gLsaUK8kiSV1nQT9lOBz2Tmf2TmFuDtFMVjLszM71U7ZeZlwOeB503jPT4IjFa2PRUOk5MkqfV2Av0RMX9M+2IeSrB3UiTnj5qgjyRJXW+6CfuvA9+ref79yuN36vT9NnB0IxuPiGcAZ1FUj53qgfsciqv+1WVac8dKktQuIuKwiPiNVscxxvUUt9w9bUz7copEHR46h1g5QR9JkrredBP2B9j3/vfqFDNjC8xUX5vynEYRcSBwHvC5zLx4qus5TE6S1IXeCNzQ6iDGuAb4GbC22hARTwSeClwOe2eWuQ54bUREpc/hwCnVPpIkaZpF54Cfsu8V7F8Crwd+WKfv0RRVX6fqLRQVYk+PiMWVtkOhqEwPPJCZQw1HLEmSZiwiTgYOGO9H9cx8ICL+GviriLgV+BbwEYpzhX+t6XoOxW1zn46I/wT+D9BD8aO9JEli+gn7Vopic0BxdRv4xDh9n0PxK/pUPY9iSHu95P8u4ArgxAa2J0mSmuedFD+kTzQK7iNAH/AOilF2Q8AZtTPGZOYXIuJo4APAq4H7gbMyc+tsBd5KI6PJlh27uHPPMEsW9LBi6SLmz5vyAERJUpeabsL+XqZwj3hE/BowCHyugW2/GXjEmLbnUVx5P5liujhJkjSLMvNm6tzSlpknTmHdUeAvIuLjwOOBb2Xm3XX6fTgi/hU4DrguM2+fadxltHnbIOs3bWdwaHhvW39vD+tWD7BqWX8LI5Mkld20EvbM/BHwoyn0uwP4vQa3/e2xbRFxZOW1rzWyLUmS1DqZeRNw0yR9bgc6MlGHIllfu3ErOaZ959AwazduZcOa5SbtkqRxTbfo3F4R8ZJmBCJJktRJRkaT9Zu275esA3vb1m/azshovR6SJDUhYQc+FxGvbcJ2xpWZ52emN3pJkqS2sWXHrn2GwY+VwODQMFt27Jq7oCRJbaUZCfv5wCcj4p31XoyIlRHxjSa8jyRJUtu4c8/4yfp0+kmSus90i87tlZmvjoi7gPdFxOLM/AuAiHg8xZQtpwH3zfR9JEmS2smSBT1N7SdJ6j4zTtgBMvPsiLgD+JuIOIJirtVXUoz2+hTwvma8jyRJUrtYsXQR/b097BwarnsfewB9vcUUb5Ik1dOMIfFV5wGXAS8DXgP8O/DEzDwrMweb+D6SJKkQ1Jl6TeUwf16wbvUAsP//pOrzdasHnI9dkjSuZlSJPzAi3kgxbctJwHcorqwfCNw60+1LkqRxfRRY2uogNL5Vy/rZsGY5fb37Dnvv6+1xSjdJ0qSaMST+BuAoYDvwqsz874h4KfAvwEUR8aLM/GUT3keSJNXIzCFgqNVxaGKrlvVz8kAfW3bs4s49wyxZUAyD98q6JGkyzUjY51MMgT8/M0cBMvMLEfFz4D+ByyLi1My8qwnvJUmS1HbmzwtWHnN4q8OQJLWZZiTsj83M/eYjycxLIuJ3gIuAK4HHN+G9JEmSJEnqCjO+h71esl7z2reAE4CHz/R9JEmSJEnqJs2sEl9XZm4Hnjnb7yNJkiRJUieZ9YQdIDOtFi9JkiRJUgPmJGGXJEkzExGLI+IvI+IDEfHMmva3RcQtEbEnIi6JiGNbGackSWqeZhSdkyRJsygi+oBrgUdWmt4eEWuAxcB64DvAz4HnAJdHxHGZeVtLgpUkSU3jFXZJksrvLUAv8BLgGcB1wPuBPwGelpkrM/MpwMnAocBftiZMSZLUTCbskiSV36nAZzLzPzJzC/B2YClwYWZ+r9opMy8DPg88rzVhSpKkZjJhlySp/H4d+F7N8+9XHr9Tp++3gaNnPSJJkjTrTNglSSq/B9i37sxw5fGXdfoOAzHrEUmSpFlnwi5JUvn9FDiy5vkvgdcDP6zT92jgrrkISpIkzS6rxEuSVH5bKYrNAZCZ9wOfGKfvcyiK0kmSpDZnwi5JUvm9l32vsNcVEb8GDAKfm/WIJEnSrDNhlySp5DLzR8CPACLiJZn5xXH63QH83lzGJkmSZo/3sEuS1F4+FxGvbXUQkiRp9pmwS5LUXs4HPhkR76z3YkSsjIhvzG1IkiRpNjgkXpKkNpKZr46Iu4D3RcTizPwLgIh4PHAOcBpwXytjlCRJzWHCLklSm8nMsyPiDuBvIuIIimneXgkk8Cngfa2MT5IkNYcJuyRJ7ek84AXAyygS9c8D787MH7c0KkmS1DTewy5JUhuJiAMj4o3ATcBJwHcoEvYDgVtbGZskSWouE3ZJktrLDcDfAncBL8zMpwIvB14IXBQRh7YyOEmS1Dwm7JIktZf5wGuA4zLzvwEy8wsUCfszgMsq97VLkqQ25z3skiS1l8dm5vDYxsy8JCJ+B7gIuBJ4/JxHJkmSmsor7JIktZF6yXrNa98CTgAePncRSZKk2WLCLklSB8nM7cAzWx2HJEmaORN2SZJKLCKe2+g6mXlrZd3faX5EkiRprpiwS5JUbpsj4tKIeEFEzJ+sc2XatxdFxBUU97NLkqQ2ZdE5SZLK7TcppnG7ELgrIr4GbKGYh30XEMAi4LEUVeKfCxwGXAI8Ze7DlSRJzWLCLklSiWXmNuB5EbESOAs4DfhDIMd0DWA38F/Ahsz8f3MaqCRJajoTdkmS2kBmXg1cXRkW/1RgADiCInG/C9gGfCczR2f6XhFxIvDXlfe4FnhFZt7UwLqXTdDlpMy8PCIeCdxe5/V/ysxXNxKvJEmdqvQJe0Q8DDgK+FlmDrU6HkmSWikzRyiGxG+Zje1HxNOAzcANwDuAFwEXRcSTM/P+KWziB8AZddp/G3glcEvl+VOBkUpb7Y8MN04zdEmSOk6pE/aIeCvwbuBQ4MGI+Cxw1hRPGCRJUuM+AvwCeHZm7oqI8yiS6DOBcydbOTPvADaObY+INwIXZOaOStPxwA8z87PNClySpE5T2irxEfFi4IPAu4DjgPcArwD+ooVhSZLUsSLiEcAJwD9n5i6AzLwH+DzFvfPT3e5qimP5+prm44Grph+tJEmdr7QJO8VB/dzM/FhmXp+ZHwAuBV7Q4rgkSepUT6Q4N7hyTPv1FAn3dL0T2JiZN9e0HQ8sj4hbIuK+iLgmIk6ZaCMRcVBELKwuwIIZxCRJUumVMmGvFNR5M/DhMS/9in3vc5MkSc1zWOXx5jHtdwGLI6Kn0Q1GxPHA04GP17Q9GlhSeb8NwFuBhwNfiYjfmmBzZwNDNcttjcYjSVI7KeU97JWCOhfXtkXEMuA5wPtbEpQkSZ0vKo/3jmkfrjwuqPnvqXo9cE1mbq1pG6G46r4hM38OEBHnA9uBNzH+UPlzKOakr1qASbukLjMymmzZsYs79wyzZEEPK5YuYv68mHxFtaVSJuy1IuJJFPet/yHwvxTFcOr1Owg4qKbJYXKSJDWmmoyPPfOrPj+IBkTEAuAlwBtq2zPzVoo6NbVteyLiS8DLxttepejs3sKzEZ6gSuoum7cNsn7TdgaHHvrttL+3h3WrB1i1rL+FkWm2lHJI/BiPAJ4MHAw8CBw4Tj+HyUmSNDM7K49HjWlfXHnc0+D2TqM4bn9piv3vBRZFxHjHeknqWpu3DbJ249Z9knWAnUPDrN24lc3bBlsUmWZT6RP2zLwyM1cAJwInA38zTtdzgN6a5cg5CVCSpM5xI3AfsHJM+3LgvswcanB7LwWuzMy7ahsj4oyI+FCd/scCuzPzgQbfp+uMjCZX33Q3X/7u7Vx9092MjGarQ5I0i0ZGk/WbtlPvX3q1bf2m7X4XdKDSD4mvyswrIuLfgNOB19Z53WFykiTNQGbeHxGXAK+IiL/LzOHK1e6XAlc0sq1KgbrfYd+p3Kp6gTdExKcz86ZK/2cBvwv854w+RBdwSKzUfbbs2LXflfVaCQwODbNlxy5WHnP43AWmWVfKK+wRsTgiPhQRS8a8tIs2+pFBkqQ29GHgMcAXImIVRQJ9FEU1dyLiyRGxJiImOx7/FtBD/QJynwV+Dnw9Is6tFJz7KsUtbe9qyqfoUA6JlbrTnXumVu9zqv3UPkqZsAO7gT8F3lFtqEz1dgrwrVYFJUlSp8vMq4A/Bp5LMWPL84H3ZuaFlS6/B1xAkYxP5CTgAeD/1XmP3cCzgeuAV1be40JgRWb+qAkfoyM5JFbqXksWTG1Wzan2U/soZcKemb+imL7tjRGxMSLWAhcBj6P+0DpJktQkmbkROBo4FXh8Zq6ree09mRmZ+ctJtvHuzHxYZt43zus3ZOapmXlIZh6RmX+YmTc29YN0mEaGxErqLCuWLqK/t2e/KTyqguLWmBVLF81lWJoDpUzYK/6aYu7WFRTD8x4OPCczv9nSqCRJ6gKZeXdmXmwSXR4OiZW61/x5wbrVA8D4826uWz3gfOwdqLQJexY+npmPy8wFmfnbmXllq+OSJElqBYfESt1t1bJ+NqxZTl/vvv/G+3p72LBmuUUnO5QF3CRJktpAdUjszqHhuvexB8WJu0Nipc61alk/Jw/0sWXHLu7cM8ySBcW/ea+sdy4TdkmSpDZQHRK7duNWAvZJ2h0SK3WP+fPCqdu6SGmHxEuSJGlfDomVpO7iFXZJkqQ24pBYSeoeJuySJEltZuyQ2F89OMr539zBLbvu5ehFB3PGykfzsAMcSClJ7c6EXZIkqY2dc9F2zvvGDkZrbmr/wEU/4DUnLOXsUwdaF5gkacZM2CVJktrUORdt51Nf37Ff+2iyt92kXZLal2OlJDXVyGhy9U138+Xv3s7VN93NyGi9yYckSTP1qwdHOe8b+yfrtc77xg5+9eDoHEUkSWo2r7BLaprN2wZZv2k7g0PDe9v6e3tYt3rAysWS1GQXXH0zk/0mOppFv1ed8BtzE5Qkqam8wi6pKTZvG2Ttxq37JOsAO4eGWbtxK5u3DbYoMknqTLfsurep/SRJ5WPCLmnGRkaT9Zu2U+9CT7Vt/abtDo+XpCY6etHBTe0nSSofE3ZJM7Zlx679rqzXSmBwaJgtO3bNXVCS1OHOWPloJpt6fV4U/SRJ7cmEXdKM3bln/GR9Ov0kSZN72AHzeM0JSyfs85oTljofuyS1MYvOSZqxJQt6mtpPkjQ11Snbxs7DPi9wHnZJ6gAm7JJmbMXSRfT39rBzaLjufewB9PX2sGLporkOTZI63tmnDvDm5z2BC66+mVt23cvRiw7mjJWP9sq6JHUAE3ZJMzZ/XrBu9QBrN24lYJ+kvXp75brVA8yf7GZLSdK0POyAeU7dJkkdyJ9eJTXFqmX9bFiznL7efYe99/X2sGHNcudhlyRJkhrkFXZJTbNqWT8nD/SxZccu7twzzJIFxTB4r6xLkiRJjTNhl9RU8+cFK485vNVhSJIkSW3PIfGSJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEkldECrA5AkSZKkTjcymmzZsYs79wyzZEEPK5YuYv68aHVYKjkTdkmSJEmaRZu3DbJ+03YGh4b3tvX39rBu9QCrlvW3MDKVnUPiJUmSJGmWbN42yNqNW/dJ1gF2Dg2zduNWNm8bbFFkagcm7JIkSZI0C0ZGk/WbtpN1Xqu2rd+0nZHRej0kE3ZJkiRJmhVbduza78p6rQQGh4bZsmPX3AWltmLCLkmSJEmz4M494yfr0+mn7mPCLkmSJEmzYMmCnqb2U/cxYZckSZKkWbBi6SL6e3sYb/K2oKgWv2LporkMS23EhF2SJO0VESdGxLURcW9EfD0ijmlw/UdGRNZZ/nFMv+Mi4vKIuCcivhMRxzf3k0hS682fF6xbPQCwX9Jefb5u9YDzsWtcJuySJAmAiHgasBk4CHgHRT2kiyLioAY281RgBPhj4IyaZW/CHhGPBi4FlgL/B7gZ2BwRi2f8ISSpZFYt62fDmuX09e477L2vt4cNa5Y7D7smdECrA6gnIuYB7wTOBPqAO4APZOaGlgYmSVJn+wjwC+DZmbkrIs4DbqQ4Hp87xW0cD/wwMz87QZ/3AAcDKzLzpoj4O+Ba4GzgzdMLXZLKa9Wyfk4e6GPLjl3cuWeYJQuKYfBeWddkynqF/T2V5WvAWcANwCcj4sUtjEmSpI4VEY8ATgD+OTN3AWTmPcDngdMa2NTxwFUTvE8Aq4EvZeZNlfcZAc5v8H0kqa3MnxesPOZwTnvKo1h5zOEm65qS0iXsEXEE8FbgTZn5isw8DziF4ir7q1oanCRJneuJFOcFV45pvx44roHtHA8sj4hbIuK+iLgmIk6peb0fWDTO+xwTEYeMt+GIOCgiFlYXYEEDcUmS1HZKl7ADhwDrgU9UGzJzmGJI3hGtCkqSpA53WOXx5jHtdwGLI2LSOYcq96YvqWxrA8UP8A8HvhIRvzWF9wE4aoK3OBsYqllumywmSZLaWekS9sy8OTPPycwHq20RcSBwLPD91kUmSVJHq47NvHdM+3DlcSpXs0coatA8LTP/KjPPBZ4F/BR4UxPe5xygt2Y5cgoxSZLUtkpZdK6OVwMLgXEL2FQq2NZWsXWYnCRJU1dNmMebeWjSSvGZeSvwwTFteyLiS8DLZvo+mXk/cP/eFcL7PyVJna10V9jHqkzx8h7gssz82gRdHSYnSdL07aw8jh2SXp1qbc8Mtn0vsKgyYm4230eSpI5S+oSd4h64QyimlJmIw+QkSZq+G4H7gJVj2pcD92Xm0GQbiIgzIuJDdV46FtidmQ9UKs//eJz3gYcSekmSul6pE/aIOBM4HXhdZt44Ud/MvD8zd1cX/IVekqQpqww3vwR4RbXAXOWK+EuBK6a4mV7gDRFxTLUhIp4F/C7w1Zp+m4AXV0bRVZ1BMX/7HdP/FJIkdZbSJuwR8XTg74HPZOb5LQ5HkqRu8GHgMcAXImIV8J8UQ9c3AETEkyNiTUSMVwPns8DPga9HxLkRcT5Foj4EvKum399T3Kv+5Yg4JSL+EXh69X0kSVKhlAl7RDyG4tf3bcCftTgcSZK6QmZeBfwx8FzgYuD5wHsz88JKl98DLgDqTvFWGeH2bOA64JWV9S8EVmTmj2r6/biyrccCFwGvAv4J+HjzP5UkSe0rMrPVMewnIq4GngG8A7i19rXM3DjFbSwEhoaGhli4cGHzg5QkqQG7d++mt7cXoLeS2JZWRBwOrABumOyWtBm+z6HAM4GfZub3prG+x3pJUqk0+3hfumndKicJz6g8/WCdLlNK2KVOMjKabNmxizv3DLNkQQ8rli5i/jynM5I0OzLzboor7LP9Pr9k33vbJUlSjdIl7JWTBDMRqWLztkHWb9rO4NDw3rb+3h7WrR5g1bL+FkYmSZIkaTaV8h52SYXN2wZZu3HrPsk6wM6hYdZu3MrmbYMtikySJEnSbDNhl0pqZDRZv2k79apMVNvWb9rOyGj56lBIkiRJmjkTdqmktuzYtd+V9VoJDA4Ns2XHrrkLSpIkSdKcKd097JIKd+4ZP1mfTj9JkiQ1hwWBNVdM2KUSqPelv2RB3WmO9zPVfpIkSZo5CwJrLpmwSy023pf+u58/QH9vDzuHhuvexx5AX2+R3EuSJGn2VQsCjz03qxYE3rBmuUm7msp72KUWmqgK/J99bisvPK74wh87wKr6fN3qAYdfSZIkzQELAqsVTNilFpnKl/6F1w3yiZf9Jn29+w577+vt8RdcSZKkOWRBYLWCQ+KlFpnql/4jDjmIK9/2HAubSJIktZAFgdUKJuxSizTypT9/XrDymMNnOSJJkiSNx4LAagWHxEtzbGQ0ufqmu7nhjl9Oqb9f+pIkSa23Yuki+nt79qstVBUUhYMtCKxm8gq7NIfqVYQfj1XgJUmSymP+vGDd6gHWbtxKwD51iCwIrNniFXZpjoxXEb4ev/QlSZLKZ9WyfjasWW5BYM0Zr7BLc2CiivD19PX2sG71gF/6kiRJJbNqWT8nD/RZEFhzwoRdmgOTVYSvet1Jx/DMxxzhl74kSVKJWRBYc8WEXZoDU60I/9hfW+CXvyRJkiTAe9ilOeE0IJIkSZIaZcIuzQGnAZEkSZLUKBN2aQ5UpwEB9kvarQgvSZIkqR4TdmmOOA2IJEmSpEZYdE6aQ04DIkmSJGmqTNilOeY0IJIkSZKmwoRdmoaR0fQquSRJkqRZZcIuNWjztkHWb9rO4NBDc6v39/awbvWA96FLkiRJahqLzkkN2LxtkLUbt+6TrAPsHBpm7catbN422KLIJEmSJHUaE3ZpikZGk/WbtpN1Xqu2rd+0nZHRej0kSZIkqTEm7NIUbdmxa78r67USGBwaZsuOXXMXlCRJkqSO5T3s0jjGFpbbOXTflNa7c8/4Sb0kSZIkTZUJu1RHvcJyiw552JTWXbKgZ7bCkiRJktRFTNilMaqF5cbeif7ze3414XoB9PUWU7xJkiRJ0kx5D7tUYyqF5eqpzsC+bvWA87FLkiRJagoTdqnGZIXlqhYdcuA+z/t6e9iwZrnzsEuSJElqGofESzWmWjDu3S94En0Le/YWpFuxdJFX1iVJkiQ1lQm7VGOqBeP6Fvaw8pjDZzkaSZIkSd3MhF2qsWLpIvp7e9g5NFz3nnULy0mS1Dxjp1B1xJok7cuEXaoxf16wbvUAazduJdi30JyF5SRJap56U6j29/awbvWANWEkqcKic9IYq5b1s2HNcvp69x0eb2E5SZKaozqF6thCrzuHhlm7cSubtw22KDJJKhevsEt1rFrWz8kDfQ7Tk9SVIuJE4K+BAeBa4BWZeVMD688D3gmcCfQBdwAfyMwNNX0eCdxeZ/V/ysxXTzt4ld5kU6gGsH7Tdk4e6PO4K6nrmbBL45g/LywsJ6nrRMTTgM3ADcA7gBcBF0XEkzPz/ilu5j0UCftngauAlwGfjIifZea/V/o8FRgBXgmM1qx744w/hEptsilUExgcGmbLjl0ehyV1vVIn7BFxBLAVOCMzL29xOJIkdYOPAL8Anp2ZuyLiPIok+kzg3MlWrhy73wq8KTM/Vmm7ALgZeBVQTdiPB36YmZ9t9gdQuU11CtWp9pOkTlbae9gj4uHAF4EjWx2LJEndICIeAZwA/HNm7gLIzHuAzwOnTXEzhwDrgU9UGzJzmCLpP6Km3/EUV9/VZaY6hepU+0lSJytlwh4Ri4H/BR7X6lgkSeoiT6Q4N7hyTPv1wHFT2UBm3pyZ52Tmg9W2iDgQOBb4fk3X44HlEXFLRNwXEddExCkTbTsiDoqIhdUFWDCVmFQu1SlUx7s7PSiqxTuFqiSVNGEHXg48CLyw1YFIktRFDqs83jym/S5gcURM95Lnq4GFFPe0ExGPBpZU3m8DxRD6hwNfiYjfmmA7ZwNDNctt04xHLVSdQhXYL2l3ClVJ2ldZE/YLgZOAu1sdiCRJXaSaId07pr16M3HDV7Qro+beA1yWmV+rNI9QFKV7Wmb+VWaeCzwL+Cnwpgk2dw7QW7N421ybcgpVSZqaUhady8wdABFT/2U1Ig4CDqppcpicJEmNqSbm4134PIjGbaC4r/3MakNm3gp8sLZTZu6JiC9RVJSvq1Klfm+l+kbOE1Q+TqEqSZMrZcI+TWcD61odhCRJbWxn5fEo4Mc17Ysrj3sa2VhEnAmcTjGP+1Sma7sXWBQRB2bmA428l9qTU6hK0sTKOiR+OhwmJ0nSzNwI3AesHNO+HLgvM4emuqGIeDrw98BnMvP8Ma+dEREfqrPascBuk3VJkgodk7Bn5v2Zubu60OBVAEmSul1lyPklwCuqBeYqFd5fClwx1e1ExGOATcA24M/qdOkF3hARx9Ss8yzgd4GvTvsDSJLUYTppSLwkSZq5DwPfAL4QERuAsyiGyL8OICKeDDwZ+Hzt1G1jXEAx5/pHgdNr7zXPzI0U1eLfAXw9Iv6Lou7Miykqv79rFj6TJEltyYRdkiTtlZlXRcQfA/9AMb3qKPDezLyw0uX3KGrGfAn45dj1I+Jw4BmVpx8c+zqwMTN3R8SzgY8Br6S4d/1C4N1TvNddkqSuYMKuUhoZTavGSlKLZObGiLgYWAHcUJtEZ+Z7KKZpG2/du9m/yny9fjcAp844WEmSOpgJu0pn87ZB1m/azuDQ8N62/t4e1q0ecF5WSZojlcT74lbHIUlSNyt10bnMvDkzIzMvb3Usmhubtw2yduPWfZJ1gJ1Dw6zduJXN2wZbFJkkSZIkza1SJ+zqLiOjyfpN28k6r1Xb1m/azshovR6SJEmS1FlM2FUaW3bs2u/Keq0EBoeG2bJj19wFJUmSJEktYsKu0rhzz/jJ+nT6SZIkSVI7M2FXaSxZ0NPUfpIkSZLUzkzYVRorli6iv7dn3LmAgqJa/Iqli+YyLEmSJElqCRN2lcb8ecG61QPA/hP4Vp+vWz3gfOySJEmSuoIJu0pl1bJ+NqxZTl/vvsPe+3p72LBmufOwS5IkSeoaB7Q6AGmsVcv6OXmgjy07dnHnnmGWLCiGwXtlXZIkSVI3MWFXKc2fF6w85vBWhyFJkiRJLeOQeEmSJEmSSsiEXZIkSZKkEjJhlyRJkiSphEzYJUmSJEkqIRN2SZIkSZJKyIRdkiRJkqQSMmGXJEmSJKmETNglSZIkSSohE3ZJkiRJkkrogFYHIEmSJLXCyGiyZccu7twzzJIFPaxYuoj586LVYUnSXibskiRJ6jqbtw2yftN2BoeG97b19/awbvUAq5b1tzAySXqIQ+IlSZLUVTZvG2Ttxq37JOsAO4eGWbtxK5u3DbYoMknal1fYu8CvHhzlgqtv5pZd93L0ooM5Y+WjedgB/lYjSZK6z8hosn7TdrLOawkEsH7Tdk4e6HN4vKSWM2HvcOdctJ3zvrGD0Zqj0gcu+gGvOWEpZ5860LrAJEmSWmDLjl37XVmvlcDg0DBbduxi5TGHz11gklSHCXuHqS2e8j/b7+Ar1+8/pGs04VNf3wFg0i5JkrrKnXvGT9Zr7dw9zNU33W1BOkktZcLeQeoVT5nIed/YwZuf9wSHx0uSpK6xZEHPlPq97yvfZ9c9D+x9bkE6Sa1gptYhxiueMpHRhAuuvnn2gpIkSSqZFUsX0d/bw2TXymuTdbAgnaTWMGHvABMVT5nMLbvubXo8kiRJZTV/XrBudXFL4NikfaIkvnqetX7TdkZGp3PWJUmNM2HvAJMVT5nI0YsObnI0kiRJ5bZqWT8b1iynr3ff4fGLDnnYhOvVFqTrZiOjydU33c2Xv3s7V990tz9gSLPIe9g7wFSLp4w1L+CMlY9ubjCSJEltYNWyfk4e6NtbrHfJgh52Dt3HX3zxuknXne65VyeoVzPJ+/ul2WPCXgK1ld2nU4V0qsVTxnrNCUstOCdJkrrW/Hmxz9RtV99095TWm+65V7ur1kwaez29en//hjXLTdqlJjNhb7Fm/EpZLZ6yc2h4Svexzwuch12SJGmMyc6pAujrLS6udJuJaiYlxb5Zv2k7Jw/0Of2d1EReXm2h8Sq7N1qFdCrFU05ffiR/tPJo3v38J/LD951isi5JkjTGVM6p1q0emDQh7cR7vCermeT9/dLs8Ap7izT7V8pq8ZSxV+v7vKdIkiRpymZ6TtWp93hP9b79br6/X5oNJuwt0sivlLX3Vk2kXvGURu+HlyRJ6nbTPafq5Hu8p3rffrfe3y/NFhP2FpmtXynHFk+RJElS4xo9p5rq6MkFBx3Iz+65v+0urHh/v9QaJuwt4q+UkiRJnWOqoydf/k/f2tvWTkPlq/f3r924lYB9kvZG7u+X1BiLzk1itoqGVH+lHO8rLSi+xP2VUpI0lyLixIi4NiLujYivR8Qxs7GNiDguIi6PiHsi4jsRcXxzPoHUGtO5d7u20HA7FKqr3t/f17vvBaW+3p62Hu4vlVlklu/LoBkiYiEwNDQ0xMKFC6e1jdkuGlK9zwnq/0rpF58kdY7du3fT29sL0JuZu1sdTz0R8TTgG8ANwD8BLwL6gCdn5v3N2kZEPBr4NvBL4O+BZwEnAE/IzJ81EO+Mj/VSs1x909384XnXNLxeAIcdfCAHHTCPnbsf+mdW5qvvI6NpzSRpHM0+3puwj2O8oiHNTqY7tZKoJGlfbZKwXw48ARjIzF0RcQhwI/DBzDy3WduIiPOBlwLLMvOmiJgPXAtcmplvbiBeE3aVxsho8qwPXTruPd6N8gKO1J6afbwv7ZD4ZgzJm67JioZAUTSkGUOVVi3r58q3PYd/e80z+NgfPIV/e80zuPJtz/GLWZI0pyLiERRXuf85M3cBZOY9wOeB05q1jYgIYDXwpcy8qdJnBDh/qu8jldFEc7hPR7PPOSW1p1Im7JXhdJuBg4B3UHxnXRQRB83F+zcy5VozVKuQnvaUR7HymMMdUiRJaoUnUpwXXDmm/XrguCZuox9YNE6fYypX5KW2NN493tPV7HNOSe2nrFXiPwL8Anh2ZTjdeRTD6c4EpjQkbyZma8o1SZJK7LDK481j2u8CFkdET2ZOduCbdBuT9AE4CvhhvY1Xfriv/fF+wSTxSHNu7Bzuiw89iDd/8bvcsfv+aQ+V95xT6l6lu8LejCF5M+WUa5KkLlQd3nXvmPZqpjCV5Hgq25jJ+5wNDNUst00hJmnO1Y6efOZjFvOeFz4JmP5Qec85pe5VuoSd5gzJmxGnXJMkdaFqwjz28Fd9PpXb0qayjZm8zzlAb81y5BRiklpuvKHy/b09HHbwgZ5zShpXGYfEH1Z5vHlM+4RD8po5TK5aNGTtxq0E9adcW7d6wHvNJUmdZGfl8SjgxzXtiyuPe5q0jZ/X9GGcPnVVpoXbO+9VUb9Oag9jh8pXp0P7n+07PeeUNK4yXmGf7lC5pg6TG++X0L7eHqfXkCR1ohuB+4CVY9qXA/dl5lAztlG5ze3H4/SBh5J+qePUKzTsOaekiZRuHvaIeC7wNeCYzPxxTfvJwCXAUZm5XzI+zhX222Y6N+vIaO73S6i/ckqSGtUm87B/ieLWtOMyczgiDgRuAr6fmac0axsR8XfAGcDjM/NnlbZrKPbNExuI13nY1TE855Q6Q7OP92UcEj+tIXmzNUyu+kuoJEld4MPAN4AvRMQG4CyK4/HrACLiycCTgc9n5oPT2UbF3wOvBr4cEe8Hfh94OvDGpn8iqU14zimpnjIOiW/GkDxJktSgzLwK+GPgucDFwPOB92bmhZUuvwdcAIxbsnoK26Aygu73gMcCFwGvAv4J+HiTP5IkSW2tdEPioWlD8hwmJ0kqjXYYEl8VEYcDK4AbMvPG2dpGRBwKPBP4aWZ+bxrv4bFeklQq3TAkHqY2nE6SJM2CzLyb4ur4rG4jM38JfHUm7yNJUicrZcKemVdFxB8D/wC8EBhlzHA6SZIkSZI6WSkTdoDM3BgRFzPDIXmSJEmSJLWj0ibs0JwheZIkSZIktaMyVomXJEmSJKnrmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEKlntatGXbv3t3qECRJ8ng0i9y3kqSyaPYxKTKzqRssi4h4FHBbq+OQJGmMIzPz9lYH0Qk81kuSSqwpx/tOTtgDeCSwp9WxNNkCipOTI+m8z9Zs7qvGuL8a4/5qjPursAD4aXbqwXeOtcGx3r/7xri/ps591Rj3V2PcX42pt7+adrzv2CHxlZ3TcVcwinMTAPZkpmMAJ+C+aoz7qzHur8a4v/bq5s/edGU/1vt33xj319S5rxrj/mqM+6sx4+yvpu03i85JkiRJklRCJuySJEmSJJWQCXv7uR9YX3nUxNxXjXF/Ncb91Rj3l7qRf/eNcX9NnfuqMe6vxri/GjOr+6tji85JkiRJktTOvMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSVkwl4SEXFiRFwbEfdGxNcj4pgG139kRGSd5R9nK+ZWiogjIuLWiDhxmuufFRE/joh7IuJfIuLg5kZYLjPZXxHxsnH+ttY0P9LWioh5EfHuyr56ICJui4i109jO6RGxvfLv+cKIOGI24m21ZuyviPitcf6+3jVbcUvNNtNjUmUbXXFcmun5TmUbHXvO06T9M+NttAvPnxvnOXRjynAOfUCjb6zmi4inAZuBG4B3AC8CLoqIJ2fmVKcHeCowArwSGK1pv7GZsZZBRDwc+CJw5DTXfz3w98BXgHOBNwCfBP6kSSGWykz3F8Xf1g+AD45p/+ZM4iqp9wDvBD4LXAW8DPhkRPwsM/99KhuIiNMo9vc3gU8Br6g8P2k2Am6x9zDD/UXx93UX8KYx7d9tUozSrGrCd2zXHJeadL4DHXrO04z908R9XHqePzfOc+jGlOYcOjNdWrwAlwM7gUWV54cAg8DrG9jGemBbqz/LHOyrxRSJwe1AAic2uP5CYIjiC35epe2pFF/ST2715yvb/qps4wrg463+LHOwr44AhoE31rT1VP5tbp7iNuYBPwauB3oqbY+sbHd1qz9j2fZXZZ1/Ab7S6s/j4jKdpUnfsV1zXGrG+U5lvY4852nS+WBT9nE7LJ4/N7y/PIeew/1V2UZTzqEdEt9iEfEI4ATgnzNzF0Bm3gN8HjitgU0dT/FH1eleDjwIvHCa659E8YXz0cwcBcjMbwPXzWCbZTaj/RUR84DldMff1iEUB+5PVBsyc5jiV/apDmk/FlgKfKKyLpn5U+BiGvv33A6asb+ge7671JlmekyCLjkuNfF8Bzrwe6MZ+6fJ+7jUPH+eFs+hG1Oac2gT9tZ7IsX/hyvHtF8PHNfAdo4HlkfELRFxX0RcExGnNCvIErmQ4gvj7mmu/ySKX8nGDkVpdH+3i5nur8cDhwKviYi7I2J3RHw5Ih7XtAhLIjNvzsxzMvPBaltEHEiRhH9/ipt5UuVxpv+eS68Z+ysiDgWeADw/InZW7oe7NCJWzE7U0tRFRG9E9E2wPJyZf8dChxyXJttfwADNOd+Bzjznacb5YLPOKduB58+N8xy6MaU5hzZhb73DKo83j2m/C1gcET2TbSAiHg0sqWxrA/BW4OHAVyLit5oUZylk5o7MHJnBJg4D7s7MX45pvws4agbbLaUm7K+nVx4DWAd8GFgJXBERh80wvHbwaopfkz87xf6HVR5vHtPekX9fdTS6v46nOA4dDHyI4m/s0cClEfEbsxGg1ICPUQyvHW95aRO+Y6FzjkuT7a/qCf3NY9ab8vkOdPQ5z2GVx5vHtDeyf5qxjXZxWOXx5jHtnj+Pw3PoxpTpHNqic60Xlcd7x7QPVx4X1Pz3eEYoCj9tyMyfA0TE+cB2ikJO3TLUZyqC/fc1FPt4wRzH0g5+RFFQ5ONZuRknIi4CrgXWAB9vYWyzKiIWUxRVuywzvzbV1SqP941p7/i/r2nur7uAtwEfy0qBoIj4N4q/u7OAv5yFUKWp+jCwcYLXpzryZjKdclyabH/9WuVxJuc70LnnPM04H2zGNtqF589zr1O+q+ZK086hTdhbr/plEmPaq88PmmwDmXkrY6oPZuaeiPgSRdVmPWSY/fc1lbZJ93W3ycyrGHPAysytEfF9imFCHZuwU/zafghwZgPrDFMMF5vHvtVmu+Hvq+H9lZnfZ0zSk5m3R8SldGZVfbWRzNxOceI+2zriuDTZ/oqI51b/c+xLlccpfdYOPueZ8flgk7bRLjx/nnsd8V01V5p5Du2Q+NbbWXkcO5RkceVxzwy2fS+wqHJfqQo7gf6ImD+mfTEz29fd5l6gr9VBzJaIOBM4HXhdZjYytctOigPXo8a0d/Tf1wz213g6+u9LGqNbjkuzeb4D7X/O04z9M9v7uEw8f5573fJdNdsaPscxYW+9GymGz64c074cuC8zhybbQEScEREfqvPSscDuzHxg5mF2jOspRpY8bUz7ch768ldFRHwkItaMaTuAolDYz1oT1eyKiKdTzDH6mcw8v8HVv1d5rPfvuSP/vmayvyLizRHx5jovHUuH/n1JdXTLcWnG5zvQ0ec8zdg/TdnHbcLz57nXLd9VTdHMc2gT9har3Ld5CfCKaoGMyi96L6WYu28qeoE3RMQx1YaIeBbwu8BXmxtx27uG4h/J2mpDRDyRYh7Jy1sUU5ktBd4eEbVDnd5CUVis4/62IuIxwCZgG/Bnja6fmT+hmN7ktRERlW0eDpxCB/59zXR/Af3AWyvT81S3+YcU1X877u9LGkdXHJeadL4DHXrO04z908R9XHqeP7dEV3xXNVHzzqFnOpG7y8wX4LcoCl98GVhFMY1AAi+svP5kiuIEB4yz/kLgp8DtwLnA+cA9wC7gca3+fLO0zx5d2Ucnjmk/GThlknXfVln3/cBq4IfAbuDIVn+usu0v4BkUc1B+G/hb4CuV7VwL9LT6c83Cfrq68vnOrvyb27tUXl8JnD7JNl5a2cZ5lX/PV1X24fJWf76y7a/K3+UvKQqz/B3F/LkPUFT9PaLVn8/FZarLeN+xldc8Lj30OSc836n06dpznsn2z2T7Zqr7uFOWme6vTv5bmmS/1f2+8ruqufuLJp5Dt3wnuOz9n7qmcuKalS+f9TWvvafSfugE6z8WuKjyRXMX8G/AY1r9uWZxf433j+dy4NpJ1p0HfJSiKFgCvwBOa/VnKvH+eh6wlaLYyM3AR4AFrf5Ms7CPDq/so7pLpc/5wM+msK23UiSeWdlvZ7b685V1f1EMrfsGxdDGnwKfBn6t1Z/PxaWRZbzv2MprHpf2/azjnu9UXu/qc56Zng9OZR930jLT/dXJf0sTfOaZnBN2zXdVk/ZXU86ho7IxlUBl6OwK4IZsTuEmTaAyBOrxwLcy8+5Wx6POEhGPoph3+LrMvL3V8Ugqv245Lnm+M7Fm7J9u2sfd9FnLolu+q8rChF2SJEmSpBKy6JwkSZIkSSVkwi5JkiRJUgmZsEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskqYkIh4fEf8aET+IiKGIuDcifhgRfxsR/a2OT5IkSeo0B7Q6AElt40igH/i/wG3Ag8CxwJnAH0TEUzLzzhbGJ0mS5kBE5GxuPzNjNrcvtZPInNV/b5I6XES8GPgi8LbM/HCr45EkSZI6hUPipS4WEQ+PiNsi4icRcdCY1/4xIkYi4g8m2cwtlcdHzE6UkiRJUncyYZe6WGbeB6wDjgLOqrZHxDnAq4DXZ+bna9eJiJ6IWBwRR0bE84BPVV66aI7CliSpq7VbXZl2i1cqE4fES10uIuYD1wFLgN8AXg18FFiXme+t0/91wLk1TTcD78rMf539aCVJUkQ8F3gncA371pV5BbAbKFVdmXaLVyoTE3ZJRMQLgE3ApcBJwMcz8w3j9D0SeAJwKPCbwAuB8zPzY3MUriRJqqPd6sq0W7xSKzgkXhKZ+RXgO8BzgC8Ab5yg722Z+bXM/FJmrgP+GPhwRJw9N9FKktR52q2uTLvFK7UrE3ZJRMRLgeMqT/dkA0NvMvN6imT/rMn6SpKk+tqtrky7xSu1K4fES12ucsDcVFkeAF4MHJuZP2hgG9cBj8nMQ2YnSkmSOl+71ZVpt3ildnRAqwOQ1DoR8XTgv4BvAi8HjgR+HzgH+P/G9O3LzJ11tnESsAy4fJbDlSSpo2XmSES8neJH9C9T1JU5t17yW/El4IfsW1dm8UTvERGHAX/eQFh/n5m7WhWv1O28wi51qYgYAL5Bcf/YiZm5u9K+AfhT4FmZ+c2a/v8X6KcoTHcL0AM8FfgD4N7KNr47l59BkqROFBFbKRLazwMvm+qtahHxZOD/Ae/JzHPG6fNoYEcD4Tw2M29sVbxStzNhl7pQRPw6xVX1+4FnZuYdNa89ErgR+E5mPrOm/SXAH1Hc634EkBSJ+/8AH8nMn8zdJ5AkqTNV6sp8jqLW1HmZeWaD618DPCozj5qN+Oq8X1vFK7UbE3ZJkiSpBNqtrky7xSu1I6vES5IkSS1Wp67Mu4BRiroyY/v2jbONal2Za2Yv0r3v1VbxSu3KK+ySJElSC7VbXZl2i1dqZybskiRJUou0W12ZdotXancm7JIkSZIklZD3sEuSJEmSVEIm7JIkSZIklZAJuyRJkiRJJWTCLkmSJElSCZmwS5IkSZJUQibskiRJkiSVkAm7JEmSJEklZMIuSZIkSVIJmbBLkiRJklRCJuySJEmSJJWQCbskSZIkSSX0/wO2Ar9gG3t2wAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1200x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 将特征的中心重置为0\n",
    "Y = data['y']\n",
    "X = data[['x3']]\n",
    "X2 = X ** 2\n",
    "# 为在Matplotlib中显示中文，设置特殊字体\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "# 正确显示负号\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "plt.rcParams.update({'font.size': 13})\n",
    "fig = plt.figure(figsize=(12, 6), dpi=100)\n",
    "ax = fig.add_subplot(1, 2, 1)\n",
    "ax.scatter(X, X2)\n",
    "ax.set_xlabel('$x3$')\n",
    "ax.set_ylabel('$x3^2$')\n",
    "ax1 = fig.add_subplot(1, 2, 2)\n",
    "center_x = X - X.mean()\n",
    "center_x2 = center_x ** 2\n",
    "ax1.scatter(center_x, center_x2)\n",
    "ax1.set_xlabel(r'$x3 - \\overline{x3}$')\n",
    "ax1.set_ylabel(r'$(x3 - \\overline{x3})^2$')\n",
    "plt.savefig('center_data.png', dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.483\n",
      "Model:                            OLS   Adj. R-squared:                  0.414\n",
      "Method:                 Least Squares   F-statistic:                     7.016\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00706\n",
      "Time:                        17:46:03   Log-Likelihood:                -23.004\n",
      "No. Observations:                  18   AIC:                             52.01\n",
      "Df Residuals:                      15   BIC:                             54.68\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.5275      0.877     -1.742      0.102      -3.397       0.342\n",
      "x3             1.1905      1.367      0.871      0.398      -1.723       4.104\n",
      "x3_squared    -0.0332      0.442     -0.075      0.941      -0.975       0.908\n",
      "==============================================================================\n",
      "Omnibus:                        4.288   Durbin-Watson:                   2.713\n",
      "Prob(Omnibus):                  0.117   Jarque-Bera (JB):                1.393\n",
      "Skew:                           0.045   Prob(JB):                        0.498\n",
      "Kurtosis:                       1.640   Cond. No.                         28.9\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "# 使用未归一化的特征建模\n",
    "X = pd.concat([X, X2], axis=1, ignore_index=True)\n",
    "X.columns = ['x3', 'x3_squared']\n",
    "X = sm.add_constant(X)\n",
    "re = train_model(X, Y)\n",
    "print(re.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.483\n",
      "Model:                            OLS   Adj. R-squared:                  0.414\n",
      "Method:                 Least Squares   F-statistic:                     7.016\n",
      "Date:                Tue, 07 Nov 2023   Prob (F-statistic):            0.00706\n",
      "Time:                        17:46:03   Log-Likelihood:                -23.004\n",
      "No. Observations:                  18   AIC:                             52.01\n",
      "Df Residuals:                      15   BIC:                             54.68\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "=====================================================================================\n",
      "                        coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------\n",
      "const                 0.0822      0.345      0.238      0.815      -0.653       0.818\n",
      "x3_center             1.0970      0.305      3.594      0.003       0.446       1.748\n",
      "x3_center_sqaured    -0.0332      0.442     -0.075      0.941      -0.975       0.908\n",
      "==============================================================================\n",
      "Omnibus:                        4.288   Durbin-Watson:                   2.713\n",
      "Prob(Omnibus):                  0.117   Jarque-Bera (JB):                1.393\n",
      "Skew:                           0.045   Prob(JB):                        0.498\n",
      "Kurtosis:                       1.640   Cond. No.                         2.89\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# 使用重置之后的特征建模\n",
    "X = pd.concat([center_x, center_x2], axis=1, ignore_index=True)\n",
    "X.columns = ['x3_center', 'x3_center_sqaured']\n",
    "X = sm.add_constant(X)\n",
    "re1 = train_model(X, Y)\n",
    "print(re1.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
